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The  revised  Common  and  FORTRAN  listings  for  the  OIL 
code  described  herein  are  as  they  existed  on  July  1,  1965. 
The  OIL  code  has  been  in  continuous  development  for  3 
years  and  in  its  presented  form  has  been  applied  success¬ 
fully  by  General  Atomic  to  the  kind  of  problems  discussed 
later  in  this  report.  However,  the  development  and 
improvement  of  the  code  are  being  continued,  so  that 
duplication  of  results  (or  even  close  agreement)  between 
problems  run  with  the  code  as  published  and  the  code  as 
it  existed  either  before  or  after  this  time  is  not  neces¬ 
sarily  to  be  expected. 

General  Atomic  has  exercised  due  care  in  preparation, 
but  does  not  warrant  the  merchantability,  accuracy,  and 
completeness  of  the  code  or  of  its  description  contained 
herein.  The  complexity  of  this  kind  of  program  precludes 
any  guarantee  to  that  effect.  Therefore,  any  user  must 
make  his  own  determination  of  the  suitability  of  the  code 
for  any  specific  use,  and  of  the  validity  of  the  information 
produced  by  use  of  the  code. 


ABSTRACT 


The  three  principal  areas  of  activity, 

1  .  Numerical  solutions  of  problems  in  impact, 

2.  Code  development  for  solving  impact  problems,  and 

3.  \nalytical  work  on  the  theory  of  the  impac*  process, 

are  reviewed,  utilizing  wherever  possible  cited  papers  which  have  been 
published  during  this  past  year  as  part  of  the  project  work.  The  investiga¬ 
tions  covered  in  these  papers  are  described  only  briefly  in  the  present  status 
report,  familiarity  with  or  availability  of  the  original  documents  being 
assumed. 

The  major  part  of  the  present  discussion  is  devoted  to  a  status 
report  of  unfinished  work  on  the  problem  of  computing  strength-dependent 
and  viscous  impact  flows.  A  computer  program  is  described  for 
generalizing  Eulerian  hydrodynamic  codes  to  include  these  effects  and 
sample  calculations  are  given. 
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I.  SUMMARY  OF  WORK  AND  INTRODUCTION 


TO  PRESENT  REPORT 


The  major  areas  of  work  in  the  present  contract  period  are  summarized 
below  to  give  the  current  status  of  each  area  and  to  cite  where  appropriate 
project  documents  that  have  been  written. 

1.1.  STUDIES  OF  IMPACT  HYDRODYNAMICS 

Work  in  this  area  is  the  subject  of  the  paper,  "On  the  Theory  of 
Hypervelocity  Impact,  "  by  J  .  M.  Walsh  and  W.  E.  Johnson,  which  appears 
in  the  Proceedings  of  the  Seventh  Hypervelocity  Impact  Symposium, 

Volume  II,  pages  1—76,  February  1965.  This  paper  is  partially  an 
exposition  of  the  thick-target  hydrodynamics  studies  given  in  our  preceding 
annual  report;  it  also  contains  a  Part  II  which  is  devoted  to  a  selection  of 
impacts  with  thin-plate  targets.  The  cited  report  is  widely  circulated, 
and  many  of  the  results  will  also  be  the  subject  of  discussion  in  a  joint 
report  under  preparation  by  General  Atomic  and  the  Ballistic  Research 
Laboratories  (see  Section  1.5).  Accordingly,  it  does  not  seem  desirable 
to  detail  the  work  on  impact  hydrodynamics  as  part  of  the  present  report. 

1.2.  CODE  DEVELOPMENT  WORK  IN  HYDRODYNAMICS 

vv 

This  effort  has  consisted  primarily  of  the  development  of  improved 
computer  programs  for  the  solution  of  hydrodynamic  flows  which  are 
functions  of  two  space  variables  and  time.  The  principal  result  is  the 
continuous  Eulerian  code,  OIL,  developed  over  the  preceding  and  present 
contract  periods  and  documented  in  the  report,  "OIL-  A  Continuous  Two- 
dimensional  Eulerian  Hydrodynamic  Code,"  by  W.  E.  Johnson,  published 
as  General  Atomic  Informal  Report  GAMD-5580  (Rev.),  January  1965. 
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Additional  improvements  that  have  been  made  in  the  OIL  code  subsequent 
to  the  time  the  above  report  was  issued  are  contained  in  the  present  report 
as  part  of  Section  III.  The  more  significant  changes  are  an  improved 
treatment  of  free  surfaces, and  a  capability  for  solving  time-dependent  x-y 
flows  in  addition  to  time -dependent  axi  symmetric  flows. 

1.3.  CODE  DEVELOPMENT  WORK  IN  VISCOUS  AND  STRENGTH- 
DEPENDENT  FLOWS 

For  most  of  the  past  year,  the  principal  objective  of  the  effort  has 
been  the  development  of  a  suitable  code  for  the  solution  of  two-dimensional 
viscous  flows  and  two-dimensional  flows  in  which  strength- dependent 
deformation  is  important.  Since  Eulerian  hydrodynamics  codes  have  been 
more  successful  t  .an  Lagrangian  codes  in  handling  the  large  material 
distortions  characteristic  of  impact,  it  was  decided  to  retain  the  Eulerian 
code,  OIL,  and  to  generalize  it  to  treat  the  tensor  forces  that  arise  from 
viscosity  and  strength. 

In  adding  viscous  and  strength  options  to  the  OIL  code,  the  approach 
used  has  been  to  leave  the  hydrodynamic  capabilities  of  the  code  essentially 
unmodified  and  to  add  a  separate  phase  in  which  viscous  and  strength  forces 
are  then  taken  into  account  to  alter  the  velocities  and  specific  internal 
energies.  In  other  words,  the  average  hydrostatic  stress  and  compressions 
are  accounted  for  in  the  hydrodynamics  as  usual;  the  task  of  the  new  phase 
(called  PH3)  is  to  take  account  of  the  stress  deviator  tensor  which  arises 
from  the  strain  de viator  tensor.  For  simplicity  and  to  circumvent  any 
need  for  storing  components  of  the  stress  or  strain  deviators,  the  consti¬ 
tutive  equations  have  so  far  been  picked  to  be  of  special  types.  Specifically, 
the  strength  is  represented  by  a  rigid-plastic  set  of  constitutive  equations 
and  the  viscosity  has  so  far  been  taken  to  be  Newtonian.  Generalization 
of  these  classical  models  will  presumably  be  straightforward,  although 
the  retention  of  acn  elastic  phase  will  require  storing  components  of  the 
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The  viscous  and  strength  generalizations  of  the  OIL  code  are  given 
in  considerable  detail  in  the  present  report.  The  basic  equations  are  the 
subject  of  discussion  as  Section  II  and  the  code  is  described  in  Section  III 
and  is  reproduced  as  Section  IV.  The  viscous  and  strength  program,  PH3, 
is  of  such  a  nature  that  it  could  be  added  to  most  multidimensional  PIC- 
type  or  Eulerian  hydrodynamic  codes  without  disturbing  the  hydrodynamics 
capability,  as  was  the  case  with  OIL.  It  is  then  possible  to  add  these 
effects  in  a  hydrodynamics  problem  or  to  omit  them  by  merely  by-passing 
PH3  each  time  step.  One  rewarding  feature  of  including  viscosity,  how¬ 
ever,  has  been  a  significantly  improved  and  smoother  hydrodynamics. 

It  seems  probable  that  some  viscosity  will  prove  desirable  in  most  purely 
hydrodynamic  problems  in  order  to  smooth  spurious  oscillations.  Examples 
of  computations  with  and  without  viscosity  are  reproduced  in  Section  III. 

While  the  viscous  option  is  apparently  giving  a  smooth  solution  in  a 
very  satisfactory  manner,  there  is  a  remaining  difficulty  in  the  form  of 
small  oscillations  in  the  strength  version  of  the  code.  These  oscillations 
(which  prevent  the  flow  from  being  completely  arrested)  can  be  ignored, 
although  additional  code  development  work  to  remove  them  will  proba  ly 
be  carried  out  prior  to  any  extensive  production  applications. 

It  is  expected  that  the  completed  code  will  provide  a  general 
capability  for  solving  strength  and  viscous  deformation  problems  and  that 
its  most  important  application  may  be  to  those  flows  where  material 
distortions  are  sufficiently  great  that  Lagrangian  schemes  become  unman¬ 
ageable.  Our  primary  attention  will  be  to  cratering  problems,  such  as 
those  occurring  in  impact,  and  some  work  of  this  type  is  presented  in 
Section  III. 

The  computing  time  for  the  strength  or  viscous  options  in  OIL  is 
roughly  equal  to  that  for  the  purely  hydrodynamic  part  of  the  calculation. 
Finally,  FORTRAN  listings  for  the  hydrodynamic  sections  of  OIL  are  also 
reproduced  in  Section  IV,  as  a  convenient  way  to  document  minor  modifi¬ 
cations  which  are  described  in  the  text. 
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I. 4.  ANALYTICAL  WORK  ON  IMPACT  AND  RELATEI  PROBLEMS 

Several  analytical  studies  pertinent  to  various  aspects  cf  impact 
mechanics  have  been  reported  during  the  past  contract  period.  These 
include : 

"Late- stage  Equivalence  and  Similarity  Theory  for  One-dimensional 
Impacts,"  by  J.  K.  Dienes.  This  paper  is  a  discussion  of  simple  ideal- 
gas  impacts  and  exposes  most  of  the  physical  principles  underlying  the  more 
complex  numerical  work  on  axisymmetric  solid- solid  impact.  Since  tne 
paper  is  presented  in  the  Proceedings  of  the  Seventh  Hypervelocity  Impact 
Symposium,  Volume  II,  pages  187—220,  February  1965,  additional  discussion 
here  does  not  appear  necessary. 

"Approximate  Treatment  of  Plane  Shock  Attenuation  in  a  Solid, "by 

J.  M.  Walsh,  General  Atomic  Informal  Report  GAMD-5214,  May  1964. 

This  is  an  approximate  theory  of  slab  impacts  of  solids  and  has  been 
previously  distributed. 

"Hydrodynamic  Flow  Equations  with  a  Plasticity  Resistance  Law," 
by  J .  K.  Dienes,  General  Atomic  Informal  Report  GAMD-5910,  December 
1964.  This  is  a  summary  of  the  relevant  theoretical  mechanics  underlying 
the  viscous  and  strength  formulations  and  is  also  largely  reproduced  as 
Section  II  of  the  present  report. 

"A  General  Form  for  Matrix  Functions  of  Nonsingular  Second-degree 
Matrices,"  by  T.  Teichmann,  General  Atomic  Report  GA-6063,  March  1965. 
This  report  does  not  apply  directly  to  impact  mechanics,  but  it  does  give 
a  theorem  in  matrix  algebra  which  was  proved  as  an  incidental  result  of 
work  by  Teichmann  on  a  similarity  theory  of  impact. 

"Impact  Crater  Size  and  Target  Strength,"  by  F.  E.  Allison,  J.  K. 
Dienes,  and  J.  M.  Walsh,  General  Atomic  Informal  Report  GAMD-6453, 

June  1965.  Suitable  simplifying  assumptions  are  used  to  show  that  the 
dependencies  of  crater  volume  on  impact  velocity  and  target  yield  strength 
are  closely  related.  For  example,  a  proportionality  of  volume  to  impact 
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velocity  leads  to  the  result  that  crater  volume  varies  inversely  as  the  first 
power  of  the  yield  strength.  The  report  has  been  distributed  and  additional 
work  is  planned. 

"Remarks  on  Similarity  Solutions  for  Hypervelocity  Impact," 

T.  Teichmann,  General  Atomic  Informal  Report  GAMD-6501,  July  1965. 
This  is  a  general  mathematical  discussion  of  the  similarity  methods  which 
have  so  far  been  used  on  the  hypervelocity  impact  problem.  The  report 
will  be  distributed  in  the  near  future  . 

1.5.  PREPARATION  OF  A  COMPREHENSIVE  EXPERIMENTAL- 

THEORETICAL  REPORT  ON  THE  DYNAMICS  OF  IMPACT 

In  conjunction  with  appropriate  members  of  the  Ballistic  Research 
Laboratories  and  the  Drexel  Institute  of  Technology,  it  has  been  decided 
that  a  joint  experimental-theoretical  report  on  impact  might  provide  a 
timely  and  integrated  picture  of  this  general  field  of  activity.  Much  of  the 
past  work  by  General  Atomic  will  be  presented  in  this  report,  which  is 
planned  for  completion  in  the  summer  of  this  year.  Accordingly,  the 
emphasis  in  the  present  report  has  been  on  those  subjects  which  will  not 
be  extensively  reviewed  in  thi~  forthcoming  discussion. 
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II.  EQUATIONS  FOR  VISCOUS  AND  STRENGTH -DEPENDENT  FLOWS 


The  flow  equations  for  the  motion  of  a  fluid  with  a  general  resistance 
law  can  be  written,  in  tensor  notation,  as 

if +  pe  =  0  • 


Du. 


-  S.. 


Dt  ij,j 


P  “  fl  +  jvl  u  )  =  (S.  .u.)  .  , 

Dt  \  2  k  k/  ij  i  ,j 


where  the  summation  convention  is  understood  and  S. .  denotes  a  general 

U 

stress  tensor,  I  is  the  internal  energy  per  unit  mass,  P  the  density,  and 
u^  the  velocity  vector,  0  =  u.  .is  the  divergence  of  the  velocity,  and 

D_  _  _9_  _9_ 

Dt  9u  ^  Ui  9x. 


denotes  the  convective  derivative.  It  is  appropriate  to  express  the  stress 

tensor,  S..,  as  the  sum  of  a  hydrodynamic  part,  -p6..,  and  a  deviator 
^  J  ^  J 


p?  rt,  a.  . . 
ij 


such  that 


S. ,  =  -pS. .  +  <r. . 
U  U  ij 


o’.  .  =  0  ,  p  =  -  —  S.  .  . 
li  3  i  l 


The  equations  of  motion  can  then  be  written  as 


DP 

5F  +  pe  *  0  ■ 


l. 


■W  Prager,  Introduction  to  Mechanics  of  Continua,  Ginn  and 
Company,  1961. 
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P 


Du. 

1 

~dT 


p  .  +  o\ .  . 

.1  iJ.J 


P 


D_ 

Dt 


+  p8 


(a.  .u.)  .  =  a-,  .e. .  +  a. .  .u.  , 
ij  i  .J  ij  iJ  iJ.J  i 


where 


e. . 
iJ 


+  u. 
J 


:> 


is  the  strain-rate  tensor. 

Now  the  average  stress,  j,  is  a  known  function  of  th  density,  P, 
and  specific  internal  energy,  I,  through  an  equation  of  ptate,  p  =  f(P,I). 

In  order  to  complete  the  description  of  the  flow,  "s  necessary  that  the 
stress  deviator  tensor  be  related  to  th;  strain  deviator  tensor  through  a 
constitutive  equation,  which  in  the  OIL  code  is  taken  to  be  of  the  general 
form 

<r. .  "  be  . €  . .  =  e. .  -  6. .  0/3  . 

iJ  ij  U  iJ  iJ 

Three  special  cases  are  of  particular  interest.  If  b  is  constant,  the  con¬ 
stitutive  equation  describes  purely  viscous  flows  and  b  is  twice  the  viscosity, 
ft,  of  the  fluid.  The  second  case, 


b  -  (2K2/E_)1/2,  E_  =  €..  t . .  , 

2  2  ij  ij 

describes  a  rigid-plastic  material  of  the  Prandtl-Reuss  type  for  which  the 
second  stress  invariant.  J_  =  a..  <r..,  can  be  shown  to  be  constant  and  equal 

2  2  ‘J  ‘J 

to  2K  ,  where  K  is  the  yield  stress  ir.  pure  shear.  The  third  case  defines 

equations  of  the  Perzyna  type  in  which  strain -rate  effects  are  accounted 

for  by  allowing  b  to  have  the  general  form  b  =  f(J^)'  ^or  which  the  function 

f  has  to  be  estimated  in  such  a  way  that  the  solution  to  the  equationc  of 

motion  agrees  with  material -test  observations.  In  the  notation  of  Perzyna, 
_ 

P.  Perzyna,  "The  Study  of  the  Dynamic  Behavior  of  Rate -s ensitive 
Plastic  Materials,  "  Division  of  Applied  Mathematics,  Brown  University, 
Technical  Report  No.  77,  May  1962. 
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b  =  nTOv^f),  wh  ere  F  =  (sfJT/K)  -  1.  In  order  to  make  use  of  the  Perzyna 

u 

equations  in  the  current  Eulerian  code,  OIL,  the  scalar  b  must  be  known 
as  a  function  of  the  second  strain  invariant,  E^,  which  in  turn  can  be 
obtained  from  the  known  velocity  field.  This  is  done  by  deriving  from  the 
constitutive  equation  the  relation 

h  =  b2(J2>  E2  ' 

which  must  be  solved  to  find  b  as  a  function  of  E^.  Using  the  Perzyna 
notation,  one  finds 

E2  =  y2  4>2(F)  . 

The  solution  of  this  equation  for  in  terms  of  E^  is  written,  symbolically, 

j2  =  k2  *2(e 2/t2)  . 


and,  in  general,  must  be  found  by  numerical  methods.  Then, 


or. . 
U 


K y 


<HE  /YZ)  . 


which  is  of  the  required  form  since  the  right-hand  side  can  be  obtained  in 
terms  of  the  velocity  field  and  appropriate  derivatives.  From  the  point  of 
view  of  the  OIL  code,  the  principal  feature  of  the  rigid -plastic  model  is 
that  the  constitutive  equation  does  not  depend  on  the  total  strain  but  only 
on  the  rate  of  strain,  which  can  be  determined  by  appropriate  operations 
on  the  velocity  field.  The  general  elastic -plastic  calculation,  which  is 
discussed  in  Ref.  3,  would  require  that  the  strain  itself  be  known  and, 
therefore,  would  entail  a  significant  incre.  se  in  memory,  in  complexity 
of  the  code,  and  in  computer  time  necessary  to  solve  impact  problems. 

3 

J.  K.  Dienes,  "Hydrodynamic  Flow  Equations  with  a  Plasticity 
Resistance  Law,"  General  Atomic.  Informal  Report  GAMD-5910, 
December  1964. 
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In  cylindrical  coordinates  the  strain -rate  tensor  is  given  by  the 
matrix 


where  u  denotes  the  radial  velocity  (u^  in  the  generalized  notation)  and  v 
denotes  the  axial  velocity  (u^  in  the  generalized  notation).  The  first  strain 
rate  invariant,  E  ,  is  also  the  diver'.once  of  the  velocity,  and  is,  in 
cylindrical  coordinates, 


E  =  6  =  e.  .=  —  +  —  +  — . 
1  ii  9r  r  8z 


The  second  strain -rate  invariant  is 


E 


2 


T  0  6.  ■) 

3  ij/ 


(  e. . 
\  iJ 


e. .  e. .  -  0^/3  , 
ij  iJ 


which,  in  cylindrical  coordinates,  becomes 


The  flow  equations  can  be  written  in  such  a  way  that  the  local  time 

derivatives  appear  on  the  left-hand  side  and  the  right-hand  side  is  the  sum 

of  three  terms.  The  first  gives  the  effect  of  convection,  the  second  the 

effect  of  hydrodynamic  forces,  and  the  third  term  gives  the  effect  of  the 

deviator  stresses.  The  first  two  terms  are  accounted  for  in  the  original 

4 

OIL  code  and  are  discussed  separately.  The  additional  increments  due 
to  the  deviator  stress  terms  are  accounted  for  in  the  strength  code  and 
4 

W.  E.  Johnson,  "OIL,  A  Continuous  Two-Dimensional  Eulerian 
Hydrodynamic  Code,  "  General  Atomic  Informal  Report  GAMD-5580, 
October  15,  1964. 
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are  the  subject  of  this  section.  The  full  equations  for  the  general  case 
are: 


ap 

at 


-u.p 


-  pe  , 


3u 

1 

P  -T—  =  -pu.u.  .  -  p  .  +  O’.  .  .  , 

at  j  i,j  ,1  ij,j 

P  (i  +  7  u.  u  )  =  -pu.  I  .  -  p0  +  (o-..u.)  .  . 

at  \  2  k  k/  j  ,j  r  1J  1  ,  J 

Expressions  for  the  co -variant  derivatives  in  general  coordinates 
are  derived  in  tensor  analysis.  The  appropriate  expressions  for  cylin- 
drical  coordinates  are  given  by  Sokolnikoff.  Denoting  by  a  6  the  increments 
due  to  stress  deviator  terms,  the  above  equations  lead  to  the  following 
expressions  for  the  effect  of  material  strength  in  cylindrical  coordinates: 

„  6u  _  l  a  a  l 

P  6t  r  frr  3r  ^rr  Qz  ffrz  r  ^00  ’ 


ov  a  x  3  1 

— —  =  - —  o-  +  - —  (j-  +  —  cr  , 

at  0r  rz  9z  zz  r  rz 


5(!  +  \  ukuJ  i 


6t 


9  9 

—  - —  [r  (uff  +  v<7  )  1  +  - —  (u o’  +  vo-  )  . 

r  3r  1  rr  rz  3z  rz  zz 


These  expressions  can  be  integrated  over  the  finite  volumes  associated 
with  each  cell  to  obtain  expressions  appropriate  for  the  effects  of  deviator 
stresses  alone.  As  a  result  of  this  calculation  for  the  element  of  volume 
in  Fig.  1,  one  finds 


6u  = 


2ir6t 

Am 


( 

r+Ar 

z+Az 

V 

(r  +  Ar  -  o-.nAzAr  / 

r<r 

Az  +  o’ 

rr 

rz 

r 

z 

\  2/  00  ( 

I.  S.  Sokolnikoff,  Mathematical  Theory  of  Elasticity  McGraw-Hill 
Book  Company,  1946. 
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Fig.  1--The  element  of  volume  appropriate  for  deriving  the 
*  flow  equations  in  cylindrical  coordinates;  the  contained  mass 

is  half  the  cell  mass,  Am 


6v  = 


Zir5t 

Am 


/ 

r+Ar 

z+Az  \ 

r<r 

Az  +  a 

(■*+)" 
z  ' 

rz 
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r 

6(I  +  FVk)  = 


2tr6t 
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r (<r  u  +  o  u  ) 
rr  r  rz  z 


r+Ar 


+  (r  +  ^-)(cr  u  +  o-  u  ) 
\  2  /  rz  r  zz  z 


z+Az 


These  ai  2  the  expressions  used  in  the  viscous  and  strength  options  of 
OIL. 


I 
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The  stresses  are  evaluated  at  the  cell  boundaries  in  the  code.  The 
flux  of  each  component  of  momentum  and  of  energy  across  each  cell 
boundary  is  added  to  one  cell  and  subtracted  from  the  other  in  such  a  way 
that  the  total  energy  and  momentum  are  conserved  in  the  finite  difference 
approximation. 


III.  OIL  WITH  STRENGTH  AND  VISCOSITY 


3.  1.  INTRODUCTION 

4 

In  the  hydrodynamic  version  of  OIL,  the  only  stress  acting  is  the 
scalar  pressure,  which  is  computed  from  a  given  function  of  density  and 
specific  internal  energy.  The  incorporation  of  material  strength  and 
viscous  forces  into  the  code  gives  rise,  however,  to  tensor  forces  which 
must  be  computed  using  the  constitutive  equations  and  then  accounted  for 
by  using  the  appropriate  momentum  and  energy  equations  discussed  in 
the  preceding  section.  These  calculations  have  been  programmed  into  a 
single  (optional)  phase  in  the  code  and  are  performed  each  time  step. 

This  phase  of  the  combined  code,  called  PH3,  is  currently  located  after 
PHI  (where  the  field  terms  in  the  hydrodynamic  equations  are  computed) 
and  prior  to  PH2  (where  material  transport  is  performed).  Provisions 
have  been  made  for  by-passing  PH3  to  perform  a  purely  hydrodynamic 
calculation  and  also  for  subcycling  PH3  in  order  to  split  the  time  step  for 
that  phase  only. 

The  additional  coding  required  for  the  present  options  and  anticipated 
coding  for  further  modifications  necessitated  a  decrease  in  the  maximum 
allowable  number  of  cells  from  3,  500  to  2,  500. 

Most  of  the  present  section  is  devoted  to  a  description  of  the  strength 
and  viscosity  programs.  Some  changes  have  been  made  to  the  basic  OIL 
code,  however,  and  these  are  also  described.  Among  these  are  options 
for  changing  physical  units  and  for  treating  flows  in  x-y  space  as  well  as 
the  axisymmetric  case.  An  u  ved  treatment  of  free- surface  motion  is 
also  incorporated. 
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3.  2.  DIFFERENCE  EQUATIONS  FOR  PH 3 

In  the  following  discussion,  we  have  assumed  that  Ax(i)  for  ail  i  is 
a  constant  andAy(j)  for  all  j  is  a  constant.  The  space  is  axisymmetric  and 
cylindrical  coordinates  are  used,  although  Ax  and  Ay  are  used  to  designate 
cell  dimensions  in  the  R  and  z  directions,  respectively,  as  depicted  in 
Fig.  2. 

z 


j 

Ay(j) 


♦  R 


KAL 

KA 

,T 

KAR 

L 

KL 

K 

KR 

*B 

KBL 

KB 

KBR 

i 

Ax(i) 
Fig.  2 


To  calculate  thf*  stresses  at  the  cell  boundaries,  we  compute  the 
velocity  gradients  at  points  B,  L,  T,  and  R,  as  follows:  The  velocity 
gradients  at  points  L  and  B  have  already  been  calculated  from  the  previous 
row  sweep  and  from  the  cell  to  left,  and  it  suffices  to  indicate  the 
procedure  at  R  and  T.  Referring  to  the  position  R, 


du  _  U(KR)  "  U(K) 
dR  Ax(i) 


=  SI  (FORTRAN  designation), 


dv 

dR 


du 

dz 


dv 

dz 


V(KR)  '  V 
Ax(i) 

u(KA)  *_  u 


V  +  V 

(KA) 


(K) 


=  S2  (FORTRAN  designation), 


(KAR)  U(KB)  *  U(KBR) 

_ 2 _ 2. 

2Ay(j) 

(KAR)  V(KB)  +  V(KBR) 

" _ 2. 

2Ay(j) 


=  S3  (FORTRAN  designation). 


=  S4  (FORTRAN  designation), 
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and 

^-=  ~.KAj  =  S10  (FORTRAN  designation). 

iA  X\l)  I  X  \  1  “  1  J 

With  these  velocity  gradients,  we  can  calculate  the  stresses  at  the  cell 
boundaries.  The  stresses  are  defined  as 


]  3  /  lj 


o\ .  =  b(  e. .  -  6 

ij  v  ij  ij 


6..  =  Kronecker  delta, 
ij 


_  dv  dv  u_ 
eaa  dR  dz  R 


For  a  rigid-plastic  material,  the  Prandtl-Reuss  equations  are  given 


by 


b  = 


2eK 


^  fcab  ^ab 


where  K  is  the  yield  strength  and 

2  2 


fcab  £ab  *  I  [(Ir)  +  (£)  +  (t)K(!M)  • 

Eight  of  the  nine  stresses  acting  on  a  cell  in  the  axisymmetric  case  are: 


T°P  a  ♦  .  Top 


Right 


Left 


Right 
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Fig.  3 
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Referring  to  Fig.  3,  we  define  the  hoop  stress,  which  is  normal  to 
the  plane  of  the  paper,  as  shown  here,  and  is  positive  in  tension,  as  follows 


where 


<r  =  b£^  =  DDVK  (FORTRAN  designation), 


tu ~  e  -  .  e  , 
33  33  3  aa 


and 


2u 


(K) 


'33  x(i)  4-  jc(i-  1)  ’ 


b  = 


r2K‘ 


A  ' 


Here  again,  K  is  the  yield  strength  and 


2  |YU(KR)  ~  U(KL)\  (W (KA)  ~  V(KB)  \  (  2u(K)  \ 

3  L\  2Ay(j)  /  \  2Ay(j)  )  Vx(i)  +  x(i-l)/ 


1  /U(KA)  "  U(KE)  ^  V(KR)"^aCLA2 


K 


2Ay(j) 


2Ax(i)  ) 


The  five  stresses  produce  forces  on  cell  K  in  the  radial  direction 
as  follows: 

R  R 

The  force  at  the  right  side  of  the  cell  is  F^  =  ir^  Ay(j)2itx(i); 


at  the  left  it  is 
at  the  top  it  is 
at  the  bottom  it  is 


*'l  =  -o-!  Ay(j)2irx(i-  1); 

*"21T  -  a21T  w(x(i)2  -  x(i-l)2); 

_  B  B  .  ...2  ..2. 

F21  =  _<r21  )• 


and  the  hoop  stress  has  a  radial  contribution  given  by 


F33  =  -cr33 


2TrAx(i)Ay(j)  . 
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From  the  equation  of  motion  then,  the  sum  of  these  five  forces  results  in 
the  change  Au  of  cell  K  as  follows: 


Au 


2frAt 


(K)  AMX(K)  j_ 


<Tj  Ay(j)x(i)  -  o- j  ^  Ay(j)x(i-1)  + --y—  (x(i)2  -x(i-l)2) 


~~  (x(i)2  -  x(i-l)2)  -  <r33Ax(i)A>(j) 


Only  four  stresses  acting  on  cell  K  produce  a  change  Av  in  the  axial 
direction: 

T  T  2  2 

The  force  at  the  top  of  the  cell  is  ir(x(i)  "  x(i-l)  ); 

B  B  2  2 

at  the  bottom  of  the  cell  it  is  F  =  -  cr  ir(x(i)  -  x(i-l)  ); 

bL  Ls  Ls 


at  the  right  it  is 


P12R  =  2Trx(i)Ay(j); 


at  the  left  it  is 


F12L  =  -<?l2L  2irx(i-l)Ay(j)  . 


These  four  forces  then  produce  a  change  of  Av  in  cell  K: 


Av 


2iTAt 


(K)  AMX(K) 


T  B 

— — — 2  —  — J  (x{i)2  -  x(i-l)2)  +  Ay(j)  (o'12R  *U)  -ir^^xd 


The  change  of  total  energy  E  of  cell  K  that  is  due  to  the  work  done 
by  these  forces  is 

M  =  S""*  FV 
dt  ^ 

over  all  four  sides  of  cell  K  or,  introducing  the  specific  internal  energy  I, 


Mlr  [I:  i  (“2  +  v2>]  =  X>v 
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Then.  sing  the  previous  values  for  the  interface  forces  and  velocities,  the 
internal  energy  I  of  cell  K  is 


AI 


At 


(K)  AMX(K) 


U(KR)  +  U(K)  ,.w  R, 

— — -z - —  {2nxu)Ay(j)o'11  } 


U(K)  +  U(KL)  ,,  1U 
- - - -  {Zitx(i- \) Ay ())(y  1 1  } 


+  - 


U(K)  +  U(KA)  .  T  ’■ 


2. 


’  {'r21  -  x(i_1)  )} 


U(K)  +  U(KB)  ,  B  ,  2  ,.  ,x2„ 

- - - {a21  it(x(i)  -  x(i-l)  )} 


.  V(KR)  +  V(K)  ,  R  A  , 

+ - - - {c  Av(i)x(i)2it} 


— ' 4V—  fr12L^yO)x(i-i)zw} 


,  V(K)  f  V(KA)  ,  T  ,  ,  2  1X2X. 

+ - z - -  (<r22  -  x(i-l) 


V(K)  +  V{KB)  ,  3  ,  ,.,2  l42„ 

- - - {o-22  ir(x(i)  -x(i-l)  )} 


A  .  A  ,  (Au(K))Z  ,(Av(k))2 

U(K)  Au(K)  V(K)  Av(K)  2.  2. 


All  the  velocities  are  those  from  the  PHI  hydrodynamics  calculation.  The 
T,  R,  B,  and  L  refer  respectively  to  the  top,  right,  bottom,  and  left 
boundary  of  the  cell  ''n  question. 
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3.3.  LOGIC  OF  PH 3 

For  either  the  strength  or  viscosity  options,  the  forces  are  determined 
from  velocity  gradients,  for  which  knowledge  of  velocities  in  neighboring 
cells  to  the  particular  one  being  treated  is  required.  Furthermore,  all  of 
the  velocities  used  in  computing  derivatives  must  correspond  to  the  same 
time.  These  requirements  necessitate  retaining  un-updated  velocities  for 
some  cells  as  well  as  the  updated  values.  These  dual  storage  requirements 
are  satisfied  without  additional  working  storage  by  equ:valencing  with 
working  storage  used  in  other  parts  of  OIL;  also,  the  sweep  through  the  grid 
is  done  by  rows  rather  than  columns  as  elsewhere  in  OIL  in  order  to  take 
advantage  of  the  fewer  number  of  cells  per  row  (52  maximum  instead  of 
100  in  a  column).  Hence,  less  storage  is  needed  to  retain  old  and  new 
velocities  for  the  two  rows  required. 

The  task  of  PH3,  in  general  terms,  is  to  compute  the  stresses  from 
the  old  velocities  and  to  use  the  stresses  in  the  difference  equations  in 
order  to  update  u,  v,  and  I.  A  special  check  feature  is  also  employed: 

If  the  sign  of  the  velocity  difference  between  two  adjacent  cells  changes 
during  a  time  step,  then  we  set  the  corresponding  driving  force  to  zero 
within  ‘he  time  step  when  this  occurs  (actually  accomplished  by  using  the 
average  driving  force  for  the  entire  step).  This  is  called  an  overshoot 
correction  and  serves  to  avoid  cell-to-cell  oscillations  in  the  velocity. 

The  need  to  save  old  and  new  velocities  for  some  cells  and  the 
provision  to  prevent  overshoot  offer  some  special  problems  in  the 
programming.  In  the  present  version,  the  procedure  which  is  used 
involves  working  on  three  rows  of  cells  in  each  sweep.  In  the  most 
advanced  row,  stresses  are  computed  from  the  old  velocities  and  tentative 
values  of  Au,  Av  are  computed.  In  the  second  row,  for  the  cell  immediately 
below,  the  tentative  Au,  Av  are  finalized  unless  overshoot  occurs.  Over¬ 
shoot  is  checked  at  the  upper  and  right  boundaries  and  for  the  hoop  stress 

Similarly,  for  the  hoop  stress,  the  driving  force  is  set  to  zero 
within  the  time  step  when  the  sense  of  the  corresponding  strain-rate  term, 
t changes. 
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and  is  prevented  by  reduction  of  the  appropi  iate  driving  stresses,  as 
described  above.  Finally,  AI  is  computed  for  this  second  row  cell,  using 
the  final  values  of  the  stresses.  This  AI  calculation  requires  old  velocities 
from  the  surrounding  cells,  necessitating  retention  of  these  velocities  for 
the  third  row  of  cells.  Special  procedures  are  used  for  those  cells  lying 
on  the  grid  boundaries. 

Notes  explaining  various  portions  of  the  code  are  also  given  in  the 
FORTRAN  listings  of  Section  IV. 

Timings  to  date  of  computations  using  the  strength  or  viscosity 
options  indicate  that  such  problems  are  a  factor  of  two  longer  in  computer 
time  than  equivalent  hydrodynamics  problems  without  these  effects. 

The  main  logic  for  the  strength  or  viscosity  it.  done  in  subroutine 
PH3.  In  addition,  PH3  refers  to  the  following  subroutines: 


GRADR 

GRADZ 

STRESR 

STRESZ 

HOOP 

DEJLTAU 

DELTAV 

ECALC 


Calculates  the  velocity  gradients  at  the  right 
boundary  of  the  cell  in  question. 

Calculates  the  velocity  gradients  at  the  top 
boundary  of  the  cell  in  question. 

Calculates  the  two  stresses  (normal  and  shear)  at 
the  right  boundary  of  the  cell  in  question. 
Calculates  the  two  stresses  (normal  and  shear)  at 
the  top  boundary  of  the  cell  in  question. 
Calculates  the  hoop  stress  for  the  cell  in 
question. 

Calculates  the  radial  acceleration  of  the  cell 
in  question  due  to  stress  forces. 

Calculates  the  axial  acceleration  of  the  cell 
in  question  due  to  the  stress  forces. 

Calculates  the  change  of  the  specific  internal 
energy  of  the  cell  in  question  due  to  the  work 
done  by  the  stress  forces. 
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In  addition  to  the  normal  input  required  for  OIL,  the  strength  version 
requires  the  following  additional  quantities: 


Location 

Symbol 

Description 

21 

AMDM 

Cutoff  for  strength  based  on  density;  if  <"AMDM  p^, 

forces  due  to  strength  are  not  applied  on  cell  K. 

25 

FeF 

Flag,  if  =  0.  ,  PH3  (the  strength  subroutine)  will  be 
called,  if  f  0,  no  strength. 

49 

i3 

The  number  of  times  to  subcycle  through  the  strength 
routine  (time  step  At/i3,  wher.-;  At  is  the  hydro  time 
step  used  in  PHI  and  PH2). 

66 

DXN 

Cutoff  for  the  stresses. 

71 

RSTOP 

The  factor  to  multiply  the  equation-of- state  constants 
by  to  convert  units  to  cgs  system. 

72 

SHELL 

The  factor  to  multiply  the  coefficient  of  pressure  in 
the  speed -of- sound  calculation. 

107 

Z(  107) 

Velocity  gradient  cutoff. 

5841 

DDXN 

Kq  =  yield  strength,  in  the  appropriate  units. 

5843 

DKE 

=  Newtonian  viscosity  coefficient,  in  the  appropriate 
units. 

8517 

TABLM 

Factor  on  dV/dz  critical,  shifts  the  point  where  the 
full  yield  strength  is  applied. 

13583 

VT 

Minimum  p  to  trigger  rezone. 

13586 

VVABOV 

Epsilon  for  energy  cutoff,  to  cut  down  on  the  pre¬ 
cursor. 

13587 

VVBLO 

Epsilon  for  velocity  cutoff,  to  cut  down  on  the  pre¬ 
cursor. 

The  code  is  currently  being  modified  to  incorporate  all  the  strength 
constants  and  flags  on  the  dump  tape.  Until  this  is  completed,  the  various 
quantities,  such  as  DDXN,  DKE,  TABLM,  VVABOV,  VVBLO,  and  VT, 
must  be  loaded  on  every  restart. 
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3.4.  OTHER  GENERALIZATIONS  OF  OIL 

Two  additional  modifications  have  been  made  to  improve  OIL.  One 
is  the  capability  of  treating  two-dimensional  x-y  flows  as  well  as  axi- 
symmetric  ones.  The  other  is  an  improved  representation  of  free  surfaces 
that  distinguishes  between  condensed  and  vaporized  materials  and  provides 
an  improved  transport  scheme  for  the  condensed  case.  These  two  modi¬ 
fications  are  described  below. 

3.4.1,  Axisymmetric  and  Plane  Flows 

Whether  a  problem  is  to  be  in  axisymmetric  or  plane  coordinates 
is  designated  by  a  single  flag.  If  CLAM  is  used  to  generate  the  problem, 
this  flag  is  the  fourth  word  of  card  number  2  (input  cards  for  CLAM)  and 
is  0.  for  axisymmetric  geometry,  and  any  number  other  than  0.  is  used 
to  designate  x-y  Cartesian  geometry. 

When  the  CLAM  code  is  used,  the  changes  required  are  listed  below. 
The  referenced  card  numbers  and  pages  are  described  in  Ref.  4. 

The  following  statements  replace  the  cards  labelled  1790  through 
1810  in  the  input  subroutine  (page  77): 

IF  (Q000FL)  3000,  3002,  3000 
3000  TAU(i)  =  Dx(  1) 

WS  =  1. 

GO  TO  1008 
3002  WSB  =  WSA 

WSA  =  x(i)  *  *  2 
TAU(i)  =  WS*  (WSA-WSB) 

1008  CONTINUE 

Replace  cards  numbered  1250  through  1260  in  the  routine  PH3 
(page  89)  by: 

IF  (Q000FL)  6000,  6001,  6000 
6000  TAM  =  WS5  *  Dy(j)/FMX 
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GO  TO  6002 

6001  TAM  =  WPIDy  *  WS5  *  Dy(j) 

6002  E  =  0.  0 

Replace  card  number  2050  in  routine  PH3  (page  91)  by: 

IF  (Q000FL)  6004,  6005,  6004 

6004  AM(n)  =  TAM  *  WSR 
GO  TO  4341 

6005  AM(N)  =  TAM  *  TX  *  WSR 

After  statement  7162  (card  number  1380)  in  routine  Output  (page  100) 
insert 

GAM  =  Q00FL 

If  the  subroutine  SETUP  is  used  to  generate  the  problem,  then  the 
flag  is  located  m  the  variable  GAM  (location  10  for  the  CARDS  routine). 
Referring  to  Eqs.  (1)  through  (4)  of  Ref.  4,  Eq.  (1)  is  changed  to 


9p  _  9pu  9pv 
9t  9r  9a 


the  two  momenta  equations  remain  unchanged;  and  for  Eq.  (4), 

9E  _  9Pu  9Pv 
P  9t  =  "  9r  9z  • 


3,4.2.  Free  Surface  Modification 

A  modification  has  been  made  to  the  scheme  of  handling  free  surfaces 
in  the  transport  routine  of  OIL,  As  reported  in  Ref.  4,  the  OIL  write-up, 
a  velocity  change  in  the  projectile  was  the  criterion  for  ensuring  that  the 
bottom  of  the  projectile  empty  properly  as  it  continues  to  move  upward. 
However,  this  scheme  was  not  operative  after  the  reflected  shock  broke 
through  the  bottom  surface  of  the  projectile. 

The  new  scheme  is,  again,  concerned  with  the  emptying  of  cells. 

The  criterion  for  applying  the  scheme  is  that  the  energy  of  the  cell  be  less 
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than  the  energy  required  to  vaporize  the  material.  Thus,  if  the  material 
is  a  gas,  no  special  modification  is  required  if  the  material  is  a  solid, 
the  transport  is  done  using  the  density  and  velocity  of  the  receptor  cell. 

Another  feature  of  the  code  is  the  ability  to  convert  to  the  cgs  units. 

Two  additional  input  numbers  are  required,  RSTOP  and  SHELL.  Their 
definitions  are  listed  in  the  revised  Common  list  in  Section  IV. 

3.  5.  TEST  PROBLEMS 

A  number  of  test  problems  have  been  computed  during  the  course  of 

the  present  code  development,  and  in  this  section  some  representative 

results  are  presented.  It  is  expected  that  extensive  application  of  the  code 

to  impact  problems,  and  the  associated  discussion  of  impact  mechanics, 

will  be  the  subject  of  future  reports. 

Two  one -dimensional  impacts  have  been  co  mputed  in  which  wax 

4 

plates  strike  wax  targets  at  a  velocity  of  3.  5  x  10  cm/sec.  One  of  these 

impacts  was  treated  using  the  unmodified  OIL  hydrodynamics  code  and  the 

5  3 

second  was  a  viscous  flow  problem  with  viscosity  =  2  X  10  ergs  sec/cm  . 
Cell  dimensions  in  the  two  problems  were  l  cm.  The  most  interesting  results 
of  the  calculations  are  the  pressure-pulse  and  velocity-pulse  profiles 
seen  in  Figs.  4  and  5.  The  viscous  version  of  the  problem  is  seen  to  be 
effective  in  eliminating  the  oscillations  which  are  characteristic  of  the 
unmodified  version.  Some  viscous  smearing  of  the  shock  front  is  also 
evident. 

Two  axisymmetric  wax-on-wax  impacts  were  also  computed,  with 

5 

an  impact  velocity  of  4  x  10  cm/sec.  One  was  run  with  the  hydrodynamic 

version  of  the  code,  which  differs  from  previous  calculations  of  this 

impact  only  in  the  improved  treatment  of  the  free  surface  (Section  3.4.) 

and  a  somewhat  coarser  zoning.  The  other  impact  was  treated  as  a 

4  3 

problem  in  viscous  flow  with  viscosity  1  x  10  ergs  sec /cm  .  Initial 

cell  size  in  both  problems  was  0.  0525  cm.  Successive  stages  of  t’.iu  viscous 


8.0 
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of  OIL.,  in  the  case  of  a  slab -geometry  impact 
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flow  are  given  in  Fig.  6,  and  the  crater  growth  in  these  two  problems  is 
compared  in  the  mass  distribution  plots  of  Figs.  7  and  8.  It  is  seen  that 
the  gross  features  of  the  mass  motion  are  in  very  good  agreement  with  the 
experimental  crater  growth  data  by  Karpov,^*  especially  for  the  viscous 
flow.  The  viscosity  problem  is  also  in  better  agreement  with  the  experi¬ 
mental  shock-pressure-attenuation  data, as  can  be  seen  from  Fig.  9.  The 
results  are  merely  illustrative  of  the  code,  however,  since  no  attempt  has 
been  made  to  formulate  a  realistic  constitutive  equation  for  v/ax. 

The  above  axisymmetric  impact  was  also  run  as  a  problem  in 
strength -dependent  deformation  using  the  Prandtl-Reuss  constitutive 
equations  described  in  previous  sections.  A  modification  to  this  repre¬ 
sentation  was  included  in  older  to  reduce  the  yield  strength  as  the  velocity 
gradients  became  srru.ll  and  thus  to  avoid  oscillations  in  the  form  of 
overcorrections.  The  results  were  very  satisfactory  at  high  pressures 
but  still  exhibited  some  oscillatory  behavior  as  pressures  became 
comparable  to  the  yield  strength,  which  was  taken  to  be  5  kbars  in  this 
exploratory  calculation. 

A  spherically  symmetric  test  problem  was  computed  in  cylindrical 
coordinates  in  order  to  test  the  strength  code  against  preferential  treat¬ 
ment  in  the  axial  and  radial  directions.  A  hot  sphere  of  plastic  (density 
3  10 

0.92  g/cm  and  specific  energy  7X  10  ergs/g,  sufficient  to  exert  a 

pressure  of  103  kbars)  was  allowed  to  explode  in  a  cold  plastic  atmosphere 

3 

(density  0.92  g/cm  ,  zero  pressure  and  energy).  A  yield  strength  of 

5  kbars  was  assumed,  and  this  yield  strength  was  again  diminished  as 

velocity  gradients  became  small.  Figures  10  and  11  bear  out  the  spherical 

character  of  the  solution.  Similar  tests  of  sphericity  have  previously  been 
4 

reported  for  the  purely  hydrodynamic  part  of  the  code. 

^B.  G.  Karpov,  "Transient  Response  of  Wax  Targets  to  Pellet 
Impact  at  4  km/sec,"  Ballistic  Research  Laboratories,  Report  No. 

1226,  October,  1963. 
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Fig.  7--Successive  stages  in  the  viscous -hydrodynamic  axisymmetric  impact 
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>--Mass  distribution  plots  for  the  purely  hydrodynamic  axisymmetric 
impact  problem.  The  dashed  profile  is  an  experimental  curve. 


F3ESSURE  (KBARS) 
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CM 


Fig.  9  - -Experimental,  hydrodynamic,  and  viscous -hydrodynamic  shock 

attenuation  curves  for  wax 


Fig.  10  - -Computed  pressure  profiles  in  the  axial  and  radial  directions,  for  the  spherically 

symmetric  problem  as  treated  in  cylindrical  coordinates 


IV.  LIST  OF  REVISED  COMMON  AND  FORTRAN  LISTINGS 


FOR  OIL  WITH  STRENGTH  AND  VISCOSITY 


In  the  following  list  of  the  revised  Common,  location  refers  to  the 
location  of  that  symbol  relative  to  the  begining  of  Common.  Since  the 
beginning  of  Common  is  assigned  the  same  location  for  each  subroutine,  a 
program  (CARDS)  is  available  for  changing  any  word  in  Common. 

If  changes  are  made  in  the  length  of  the  dimensional  arrays,  it  will 
be  necessary  to  change  the  locations  in  the  following  Common  list. 

The  revised  FORTRAN  listing  for  OIL  follows  the  revised  Common. 
Note  that  the  C' s  denoting  explanatory  remarks  are  at  the  far  left  margin. 
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REVTSED  COMMON 


Symbol 

Location 

No.  of 
Words 

Units 

Description 

XX 

151 

53 

cm 

XX(2)  =  X(l) 

UR 

205 

200 

None 

Note  the  many  equivalence  statements 

PR 

405 

200 

None 

Note  the  many  equivalence  statements 

YY 

605 

101 

cm 

YY(2)  =  Y(l) 

AID 

706 

i 

None 

Not  used,  this  is  a  single  material 
code 

AIX 

707 

2500 

jerks  or 
ergs/gram 

Specific  internal  energy  (x)  for  cell 
(K) 

AM 

3207 

130 

None 

Mass  of  particle  (N)  for  SHELL  code 

AMD 

3337 

1 

None 

Not  used,  this  is  a  single  material 
code 

AMX 

3338 

2500 

grams 

Total  (X)  mass  in  cell  (K) 

AREA  5838  1  None  Tag,  used  in  FH2 

BIG  5839  1  sh^or  sec”1  =  dV/dZ  critical,  computed  in  PH3 

BOUNCE  58UO  1  None  Tag  used  in  FH2 

DDXN  584 1  1  ergs  or  ^  =  yield  strength  for  the  material 

jerks/cnr 

DDVK  5842  1  ergs  or  ^  *  hoop  stress  for  cell  (K) 

jerks /cm-3 

DKE  5843  1  ergs  or  jerks  =  T)^  =  the  coefficient  of  viscosity 

/cm-'  sec  or  sh 

DVK  5844  1  sec"2  or  sh”2  =  DDXN  /[pQ  DX2] 


DX 

5845 

52 

cm 

DX(i)  =  X(i)  -  X(i-l) 

DY 

5897 

100 

cm 

DY(j)  =  Y(j)  -  Y(j-l) 

E 

5997 

1 

None 

Used  in  hoop  routine 

FD 

5998 

1 

None 

Used  in  hoop  routine 

FS 

5999 

1 

None 

Flag  in  PH2 

FX 

6000 

1 

None 

Not  used 

OUT 

6001 

1 

None 

Tag  in  particle  FH2 

P 

6002 

2500 

jerks  or 
ergs/cm3 

Material  pressure  in  cell  (K) 

PABOVE 

8502 

1 

jerks  or 
ergs/cm3 

=  CP(K)  +  P(cell  alove)]/2 

PBLO 

8503 

1 

jerks  or 
ergs/cnP 

=  [P(K)  +  P(cell  below) ]/2 

PIDTS 

8504 

1 

l/cm  sh 

l,/[At  *  DY(j)]  in  PHI, 
l./*At  in  PH2 

PPABOV  8505  1  None  Not  used 
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PRR 

8506 

1 

jerks  or 
ergs/cnw 

RJL 

8507 

1 

None 

QPT 

8508 

1 

None 

RC 

8509 

1 

cm 

ReZ 

8510 

1 

None 

RHO 

8511 

1 

gm/cm^ 

RL 

8512 

1 

None 

RR 

8513 

1 

cm 

FIG 

85lU 

1 

cm 

QOOOFL 

8515 

1 

None 

SWITCH 

8516 

1 

None 

TABIM 

8517 

1 

None 

TAU 

8518 

52 

2 

cm 

TAUDTS 

8570 

1 

2 

cm  (sh  or  sec 

TAUDTX 

8571 

1 

None 

U 

8572 

2500 

cm/sh  or/ 
sec 

UK 

11072 

1 

cm/sh  or/ 
sec 

URR 

11073 

1 

cm2/sh  or/ 
sec 

UT 

11074 

1 

None 

UU 

11075 

1 

sh  or  sec 

UUU 

11076 

1 

None 

UTEF 

11077 

1 

cm/sh  or/ 
sec 

UVMAX 

11078 

1 

, -1 

sh  or 
sec"-'- 

V 

11079 

2500 

cm/sh  or/ 
sec 

VABOVE 

13579 

1 

cm/sh  or/ 
sec 

VBLO 

13580 

1 

Cui/sh  or/ 
sec 

VEL 

13581 

1 

None 

=  [ P(K)  +  P(cell  to  the  right ) 3/2 

Not  used 
Not  used 

[X(i)  +  X(i-l)]/2.  in  PHI 

If  material  leaves  grid  in  PH2, 

ReZ  set  *  1. 

Density  of  material  in  a  cell 
Not  used 

=  [X(i)  +  X(l+l)]/2.  in  PHI 
Minimum  AX  or  AY  in  CDT  routine 
Not  used 
Not  used 

Factor  on  dV/dZ  critical 

=  n  (X(i)2  -  X(i-l)2)  =  area  in  Z 
direction 

=  TAU(i)  At  in  PHI 

Not  used 

=  R  component  of  velocity  in  cell  (K) 

=  R  component  of  velocity  in  cell  (K) 
used  in  SHELL  transport 

=  [U(K)  RC  +  U(K+1)  RRj/2. 

Signal  in  PHI,  decrease  At  next  pass 

New  At  in  PHI,  for  integrating 
backwards 

Used  in  PH3 

R  velocity  component  used  to  move 
particles  in  SHELL 

| Max  velocity  j  /Min  (Ax  or  AY) 

Axial  (Z)  component  of  velocity  for 
cell  (K) 

[V(K)  +  V(cell  above) ]/2. 

[V(K)  +  V(cell  below )]/2. 

Used  as  tag  in  PHI  and  PH2 
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VK 

13582 

1 

cm/sh  or./ 
sec 

VT 

13583 

1 

gram/cm 

VTEF 

13584 

1 

cm/sh  or 
sec 

W 

13585 

1 

None 

WABOV 

13586 

1 

jerks /gram  or 
ergs /gram 

WBLO 

13587 

1 

cm/sh  or 
sec 

W2 

13588 

1 

None 

W3 

13589 

1 

None 

WPS 

13590 

1 

ws 

13591 

1 

WSA 

13592 

1 

WSB 

13593 

1 

WSC 

13594 

1 

XL 

13595 

130 

cm 

XLF 

13725 

1 

None 

XN 

13726 

1 

cm 

XR 

13^27 

1 

None 

YL 

13728 

130 

cm 

YLW 

13858 

1 

None 

YN 

13859 

1 

cm 

YU 

13860 

1 

None 

ZMAX 

13861 

1 

cm 

i 

13862 

1 

ii 

13863 

1 

iN 

13864 

1 

iR 

13865 

1 

iWS 

13866 

1 

Axial  component  of  velocity  in  cell 
(K)  for  SHELL 

Rezone  trigger  set  if  mass  of 
P=VT  leaves  the  grid 

Z  velocity  component  used  to  move 
particles  in  SHELL 

Used  in  FH3 

Minimum  specific  internal  energy 
allowed  as  a  result  of  FH3 

Minimum  velocity  allowed  as  a  result 
of  FH3 

Not  used 

Not  used 

Working  Storage 


R  coordinate  of  particle  N 

Used  in  velocity  weighting  for 
”1  transport 

R  coordinate  of  particle  N  at  cycle 
(n-1) 

Used  in  velocity  weighting  for  SHELL 
transport 

Z  coordinate  of  particle  N 

Used  in  velocity  weighting  for  SHELL 
transport 

Z  coordinate  of  particle  N  at  cycle 
(n-l) 

Used  in  velocity  weighting  for  SHELL 
transport 

The  largest  Z  value  of  any  particle 
in  SHELL  transport 

Indices  (working  storage) 
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iWSA 

13867 

1 

iWSB 

13868 

1 

iWSC 

13869 

1 

iwi 

13870 

130 

None 

3 

14000 

1 

JN 

14001 

1 

JP 

14002 

1 

4H 

14003 

1 

K 

14004 

1 

None 

KN 

14005 

1 

KP 

14006 

1 

KR 

14007 

1 

KRM 

14008 

1 

L 

14009 

1 

M 

i4oio 

1 

MA 

i4oil 

1 

MB 

14012 

1 

MC 

14013 

1 

MD 

l4oi4 

1 

ME 

14015 

1 

MZ 

l4oi6 

1 

None 

N 

14017 

1 

NK 

14018 

1 

NKMAX 

14019 

1 

NK1 

14020 

1 

NO 

14021 

1 

NR 

14022 

1 

None 

iW2 

14023 

130 

None 

X 

152 

53 

cm 

Indices  (working  storage) 

4/ 

i  of  the  cell  (K)  where  particle  (N) 
is  for  SHELL  transport 

Indices  (working  storage) 


* 

Index  of  cell  defined  such  that 
K  -  (j-l)  iMAX  +  i+1 

Indices  (working  storage) 
i 


* 

Set  by  input,  for  length  of  Z  block, 
also  used  in  EDIT 

Indices  (working  storage) 


* 

Maximum  number  of  radiation  cycles/ 
hydro  NR  *  NRM 

=  (j)  value  of  cell  (K)  where  particle 
(N)  is  used  in  SHELL  transport 

X(i)  =  right  dimension  of  zone  (i,j) 
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UL 

205 

200 

FLEFT 

205 

200 

YAMC 

304 

100 

SIGC 

504 

100 

PL 

405 

200 

GAMC 

405 

100 

TAB 

205 

15 

AMK 

220 

15 

PK 

235 

15 

QK 

250 

15 

Y 

606 

100 

cm 

ASN 

3207 

52 

AST 

3259 

52 

ASNB 

13595 

52 

Note 

ASTB 

13647 

52 

the 

RSN 

13728 

52 

■  equivalence 

RST 

13780 

52 

statements 

SIG33 

13870 

52 

DUDOT 

13922 

52 

DVDOT 

14023 

52 

DAIX 

14075 

52  . 

KL 

14009 

1 

KAR 

l4oio 

1 

KA 

14011 

1 

KAL 

14012 

1 

KBR 

14013 

1 

KB 

l4oi4 

1 

KBL 

14015 

1 

Note  the  equivalence  statements 


Y(j )  =  top  dimension  of  zone  (i,j) 

Array  for  normal  stress  at  top 
Array  for  shear  stress  at  top 
Array  for  normal  stress  at  bottom 
Array  for  shear  stress  at  bottom 
Array  for  normal  stress  at  the  right 
Array  for  shear  stress  at  the  right 
Array  for  hoop  stress 

Array  for  the  R  component  of  acceleration 
Array  for  the  Z  component  of  acceleration 
Not  used 

Note  the  equivalence  statements 
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Location 

Symbol 

Units 

Description 

z(D 

PROB 

- 

Problem  number  (if  positive,  this  is  an  OIL  run; 

If  negative,  this  is  a  SHELL  run.) 

2(2) 

CYCLE 

- 

Cycle  number  (floating  point  value) 

z(j) 

DT 

sh.  or 

At  hydro  =  tr  -  t^n_1^ 

sec 

z(4) 

PRINTS 

- 

Cycle  frequency  for  short  print 

Z(5) 

FRINTL 

- 

Cycle  frequency  for  long  print 

z(6) 

DUMPT7 

- 

Cycle  frequency  for  binary  tape  dumps 

Z(7) 

CSTOP 

- 

Cycle  number  at  which  problem  will  stop 

Z(8) 

PIDY 

- 

it  =  3-1415927 

z(9) 

TMZ 

gm 

Total  (x  +  •)  mass  at  t  =  0  (calculated  in 

CLAM  code.) 

2(10) 

GAM 

- 

If  =  0.  (cylindrical  geometry);  otherwise 

Cartesian 

Z(ll) 

GAMD 

- 

l./(y.-l.)  computed  in  input  routine 

Z(12) 

GAMX 

- 

l./(y  -  1.)  computed  in  input  routine 

z(l3) 

ETH 

jerk  or 

Total  energy  (computed  in  CLAM  for  t  =  0). 

erg 

Changed  in  PHI  at  transmittive  boundaries  and 
in  PH2  as  mass  leaves  the  grid. 

n 

Z(l4) 

FFA 

- 

Upper  limit  for  stability  and  used  to  calculate 

n 

At,  only  if  CABLN  -  0. 

z(i5) 

FFB 

- 

Lower  limit  for  stability  and  used  to  calculate 

At,  only  if  CABIN  =  0. 

z(i6) 

TMDZ 

gm 

Total  (•)  mass  at  t  =  0,  calculated  in  CLAM  code 

Z(17) 

TMXZ 

gm 

Total  (x)  mass  at  t  =  0,  calculated  in  CLAM  code 

Z(18) 

XMAX 

cm 

=  x(iMAX) 

z(i9) 

TXMAX 

cm 

2(XMAX)  at  t  =  0,  calculated  in  CLAM  code 

Z(20) 

TYMAX 

cm 

2(YMAX)  at  t  =  0,  calculated  in  CLAM  code 

Z(21) 

AMDM 

If  the  density  of  a  cell  is  less  than  AMDM  times 
the  initial  density  (pn)>  strength  is  bypassed 
for  this  cell 

CVl 

CVJ 

' - ' 

t-3 

AMXM 

gm 

Minimum  particle  (x)  mass/2,  calculated  in  CLAM 

Z(23) 

DNN 

- 

(ETH  -  E)N"NPC/ETH 

z(24) 

DMIN 

- 

If  (ECK)  >  DMIN,  problem  will  stop  and  the  edit 
routine  will  call  dump 

z(25) 

FeF 

- 

Flag  for  omitting  strength 

z(26) 

DTNA 

sh.  or 

Atn-1 

sec 

4! 


Location 

Symbol 

Units 

Description 

Z(27) 

CVIS 

- 

If  <  0,  bottom  boundary  is  transmittive; 
otherwise  it  is  reflective. 

Z(28) 

NPR 

- 

Index  (working  storage) 

Z(29) 

NPRi 

- 

Index  (working  storage) 

Z(30) 

NC 

- 

Cycle  number  (fixed  point) 

Z(3l) 

NPC 

- 

Number  of  cycles  between  short  prints 

Z(32) 

NRC 

- 

Index 

z(33) 

iMAX 

- 

Maximum  number  of  zones  in  the  R  direction 

z(34) 

1MAXA 

- 

iMAX  +  1 

z(3b) 

JMAX 

- 

Maximum  number  of  zones  in  the  Z  direction 

z(36) 

JMAXA 

- 

jMAX  +  1 

Z(37) 

KMAX 

- 

(iMAX)(  jMAX)  +  1 

Z(38) 

KMAXA 

- 

KMAX  +  1 

Z(39) 

IMAX 

- 

Total  number  of  particles  +  1,  generated  in 
CLAM,  for  SHELL  problems  only. 

Z(4o) 

ND 

- 

Total  number  of  (•)  particles  +  1,  generated 
in  CLAM 

z(  4i) 

KDT 

- 

If  =  0,  At  has  changed,  if  /  0,  At  remains 
constant 

Z(42) 

ixMAX 

r 

Not  used 

Z(43) 

NOD 

Index 

z(44) 

NOPP 

- 

Index 

z(45) 

NiMAX 

- 

New  iMAX  before  adding  new  zones 

Z  (46) 

NjMAX 

- 

New  JMAX  before  adding  new  zones 

Z(47) 

il 

- 

Maximum  i  disturbance 

Z(48) 

12 

- 

Maximum  j  disturbance 

Z(49) 

13 

- 

Not  used 

Z(50) 

14 

- 

Not  used 

z(5D 

Nl 

t 

- 

Scratch  tape  number  for  particles  if  this  is  a 
SHELL  run 

z(52) 

N2 

- 

Scratch  tape  number  for  particles  if  this  is  a 
SHELL  run 

Z(53) 

N3 

- 

Number  of  particle  records  generated  if  this  is 
a  SHELL  run 

z(5M 

n4 

- 

Number  of  particles  -  1  per  record  (MAX  =  127) 
if  this  is  a  SHELL  run 

z(55) 

N5 

- 

Not  used 
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Location 

Symbol 

Units 

Description 

Z(56) 

lib 

- 

Number  of  particles  on  last  particle  record 
if  this  is  a  SHELL  run 

z(57) 

N7 

- 

Not  used 

Z(58) 

N8 

- 

Not  used 

z(59) 

N9 

- 

Not  used 

z(6o) 

N10 

- 

=  i  value  of  zone  that  is  controlling  At 

z(6i) 

Nil 

- 

=  j  value  of  zone  that  is  controlling  At 

Z(62) 

NRM 

- 

=  maximum  number  of  radiation  cycles/hydro 
cycle,  input  number 

z(63. 

TRAD 

sh. 

=  NR- At  RAD  =  At  HYDRO 

2(64) 

XNRG 

jerk 
or  erg 

Tota?.  energy  of  (x)  material 

z(65) 

SN 

- 

If  =  0,  code  will  decrease  At  to  correct  for 

I  <  0,  otherwise  those  I  <  0  are  left  alone 

Z(66) 

DXN 

- 

Cutoff  for  the  stresses 

Z(67) 

RADER 

gm- 
cm/  sh 

Total  positive  radial  momentum  ((x)  only) 

z(68) 

RADET 

11 

Total  positive  axial  momentum  ((x)  only) 

2, (6s) 

RADEB 

IJ 

Total  positive  radial  momentum  (x)  for  material 
under  the  target 

z(70) 

DTRAD 

- 

Not  used 

Z(7l) 

REZFCT 

- 

If  =  0,  1H2  will  not  trigger  rezone. 

z(72) 

RSTOP 

- 

Factor  for  converting  uni  „s  for  .gy 

Z(73) 

SHELL 

- 

Factor  for  converting  units  for  speed  of  ^.ound 
calculation 

Z(74) 

BBOUNL 

- 

Not  used  in  this  version 

Z(75) 

T OZONE 

gm/cm^ 

Minimum  density  for  mass  flow  at  the  free 
surface.  The  mass  flux  is  held  up  unless  it 
produces  a  density  that  is  >  than  TOZONE 

z(76) 

ECK 

- 

i  (ETH  -  E\N  /ETH  -  ExN-NFCl  / 

ETH  '  ETH  '  J/NPC 

2(77) 

SBOUND 

- 

Fraction  of  A  in  mass  weighting  velocity 

2(73) 

XI 

cm/  sh^ 

Acceleration  (R  direction)  due  to  the  stresses 

z(  79) 

X2 

cm/sh 

Acceleration  (Z  direction)  due  to  the  stresses 

Z(80) 

>Y1 

- 

Not  used 

Z(81) 

Y2 

) 

I 

Not  used 
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Location  Symbol  Units 
Z(8S)  C4BLN 

Caution:  You  must  load  a 
At  for  this  option. 

This  holds  if  SN  /  0 


z(83) 

vise 

jk/g 

z(84) 

T 

sh  or 

sec 

2(85) 

GMAX 

- 

Z(86) 

WSGD 

- 

Z(87) 

WSGX 

- 

z(83) 

GMADR 

- 

Z(89) 

GMAXR 

- 

Z(90) 

31 

, -1 

sh  or 
sec-1 

z(9i) 

S2 

If 

Z(92) 

S3 

If 

Z(93) 

s4 

If 

z(94) 

S5 

If 

2(95) 

s6 

If 

z(96) 

ST 

If 

2(97) 

s8 

II 

Z(98) 

s$ 

If 

Z(99) 

S10 

If 

Z(100) 

g 

Z(101) 

jk  or 
erg 

Z( 102 ) 

g-cm/sh 

Z(103) 

II 

Z(104) 

jk  or 
erg 

z(i05) 

SKL 

jk  or 
erg/ cm 

z(io6) 

STL 

m 

z(io7) 

. -1 

nh  _^or 
sec 

Description 


If  <  0,  code  controls  At  but,  at  Z(l39)  of 
instability 

If  =  0,  code  controls  the  At,  decreasing  At  if 

or  exceed  FFA  and  increasing  At  if 
Axl  1  ^  1  less  than  FFB 

If  greater  than  0,  At  will  remain  constant 

The  change  (pi''  ^ue  to  the  stresses 

Total  time  up  uo  cycle  N,  tn  =  tn  ^  +  At 


Maximum  of  y  or  y 

X  < 


v  and  (v  -1)  in  the  CDT  routine 
Tx  ' Tmax 

Y  /(v  _1*) 

•  i  i 

yx  '(y.  -i.) 

du/dR  \ 

]  Velocity  gradients 

dv/dR  I  the  right  boundary 
du/dZ  j  Qf  cell 

dv/dZ 

u/R  J 

du/dZ  \ 
dv/dZ  I 

du/dR  (  Velocity  gradients  at  the  top 

dv/dR  J  bouMary  of  cel1  (K) 

u/R 

Mass  thrown  away  (PH2)  if  any  cell  has  a 
p  <  T0Z0NE 

Total  energy  of  this  mass  thrown  away. 


Total  radial  momentum  of  this  mass  thrown  away 

Total  axial  momentum  of  this  mass  thrown  away 

Energy  (internal)  added  to  system  when  the 
internal  is  set  to  0.  if  I  <  0.  (PH2) 

The  normal  stress  at  the  left  boundary  of 
cell  (K) 

The  s.iear  stress  at  the  left  boundary  of  cell  (K) 
Vel.cr  .ty  gradient  cutoff 
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Location 

Symbol 

Units 

Description 

Z(lO-3) 

- 

Not  used 

z(ioo) 

- 

Not  used 

Z(ilO) 

jk  or 

Critical  en^  gy  same  value  as  Z(122)  used 

erg/g 

in  PH2 

Z(lll) 

g/cm3 

Initial  density  (pQ)  of  the  material 

Z(1I2) 

cm/sh 

Initial  velocity  of  the  projectile 

or  sec 

Z(113) 

- 

Not  used 

z(il4) 

- 

Not  used 

z(ii5) 

g/em3 

Density  (pQ) 

z(il6) 

- 

a  \ 

z(HT) 

jk/e 

Eo 

z(ll8) 

- 

b 

Z(119) 

jk/cm3 

A 

Z(l20) 

- 

v0 

L 

Z(12l) 

o 

»  For  equation  of  state 

- 

- 

Z(122) 

jk/g 

Es 

Z(123) 

- 

a 

Z(l24) 

- 

8  i 

Z(125) 

- 

- 

Z(126) 

jk/cn3 

B  / 

Z(127) 

SSI 

Z(.’°8) 

SS2 

Z(129) 

SS3 

- 

z(130) 

SS4 

- 

1 

z(i3l) 

SS5 

Z(132) 

SS6 

- 

S  Not  used 

Z(133) 

SS7 

— 

• 

z(134) 

ss8 

- 

Z(135) 

SS9 

z(136) 

SSIO 

- 

z(137) 

SS11 

Z(138) 

g/cm3 

Density  check;  if  p(K)  <  Z(133),  the  stabilit 
check  for  cell  (K)  is  bypassed 
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Location 

Symbol 

Units 

Description 

Z(139) 

- 

Percent  of  instability,  used  in  CDT  if 

CABLN  <0 

z(i4o) 

sm 

jk  cr  , 
erg, 

The  normal  stress  at  the  right  boundary  of 
cell  K 

Z(l4l) 

3TR 

ii 

The  shear  stress  at  the  right  boundary  of 
cell  K 

Z(l42) 

SNT 

ir 

The  normal  stress  at  the  top  boundary  of 
cell  K 

z(l43) 

STT 

ii 

The  shear  stress  at  the  top  boundary  of  cell  K 

z(l44) 

SNB 

ii 

The  normal  stress  at  the  bottom  of  cell  K 

Z(l45) 

STB 

ii 

The  shear  stress  at  the  bottom  of  cell  K 

Z(l46) 

- 

Not  used 

Z(l47) 

- 

The  j  interface  (projectile-target) 

z(l4C>) 

A 

10  ^cm/  sec 

Z(l49) 

B 

- 

C  =  A  +  BPe  where  A  =  and  P  is 

pressure  in  megabars 

Z(l50) 

e 

- 

REVISED  FORTRAN  LISTING 


FOR  OIL 


c-.  <r>  o  o  '  o  r» 
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IBFTC  MAIN  LIST .UfcCK »REF 


U  i  M 

E  N  S 

1  0 

N 

PH* 

00  2 1 

Ph2 

003i 

OiMfcNSiUN  AM  i  130) » 

XL ( 130) « 

YL i 130) « 

PH2 

COAL 

1  0(2300), V(*5C0),AMX 

(2500) >A 1 X ( *300 ) , 

*  Pl*50G)  t 

3  1K1(  130  *  1K2U30), 

PH2 

C07L 

4DX(32J,  X(53), 

XX(34) , 

DY(IGO), 

Y i 100 ) « 

YY  101), 

PH2 

OQSl 

5TAB( 13)  «  AMKU3), 

PK(15). 

0K(1 3). 

21130) t 

J  .ISO), 

PH2 

0090 

6T  AU( 52) f  PL (200) , 

PR (200 ) * 

UL(20C) # 

OKI 200  , 

PH2 

0100 

7FLEFT  ( ICO)  »  YAMC(IOO) 

*  SIGH  IOC)  • 

GAMC ( ICO) 

PH2 

Clio 

UIMENSION  ASN  ( 52)  «AST(32),ASNB(32)  , AST 8 i 52)  «RSN(Si), 


1RSI(52)»SIGJ3(52),UUU0T(52)  ,DVULI  <  5*)  .DAI  XI 3*:) 


COMMON 

Z  ,XX 

,  UR 

,PR  ,  YY 

PH* 

C12C 

COMMON 

AID  , AiX 

,  AM 

•AMU  , AMX 

•  AREA 

PH* 

CliC 

COMMON 

BIG  .BOUNCE 

,UOXN 

, UUVK  , UKt 

,OVK 

PH2 

Cl  AL 

COMMON 

OX  ,CY 

,FG  ,FS 

•  Fx 

PH2 

CISC 

COMMON 

OUT  ,P 

•PAbOVE 

,PBLO  ,P 1UTS  » PPAbGV 

PH2 

Cl6C 

COMMON 

PRR  , PUL 

•  CUT 

,RC  ,Kt2 

,  KHG 

PH2 

0170 

COMMON 

RL ,RR ,S I6.G0G0FL.SK1TCH 

, TABLM 

PH2 

CISC 

COMMON 

TAU  ,  TAUOT  S 

•7AUUTX 

,U  ,UK 

•  URR 

PH2 

Cl9C 

COMMON 

UT  *UU 

•ouu 

»UTEF  , U  VMAX  ,V 

PH2 

02CC 

COMMON 

V ABO Vt  , VBLU 

•  Vfc^. 

,VK  , V I 

•  VTEF 

PH2 

0210 

COMMON 

VV  , VVABCV 

•  WBLO 

«K*  ,R3 

•  KPS 

PH2 

C22C 

COMMON 

US  * KSA 

,fcS8 

,MSC  ,XL 

,XLf 

PH2 

0230 

COMMON 

XN  ,  XR 

*  YL 

,  Y  L  w  ,  Y  N 

•  YU 

PH2 

0240 

COMMON 

2MAX  ,1 

til 

,1N  ,IR 

•  IKS 

PH2 

0250 

COMMON 

iWSA  tiWSB 

•  IKSC 

•1W1  ,J 

•  JN 

PH2 

0260 

COMMON 

JP  ,  JR 

,K 

,KN  ,KP 

»kr 

PH2 

02  70 

COMMON 

KRM  ,L 

*M 

«MA  ,  ME 

•  MC 

PH2 

0280 

COMMON 

MU  ,  MB 

•  M2 

,N  ,NK 

•NKMAX 

PE2 

0290 

COMMON 

NK1  , NO 

•NR 

,iK2 

PH2 

03C0 

PH2 

0390 

PH2 

0400 

E 

Q  u  i  v 

A  L 

EMC 

E 

PH2 

0410 

PH2 

042C 

OEUUI VALENCE 

(2,12, PROB) , 

(2(2) .CYCLE), 

(2(3 ) .01 ) , 

PH2 

0430 

HZ  1  A)  ,  PRINT  Si, 

(2(5) ,PR1NTL) 

,  (2(6) 

•0UMPT7), 

(2(7) tCSTOP) , 

PH2 

0440 

212(B) ,PIUY) a 

(2(9), TMZ), 

(Zi 10) , GAM) , 

(2(11) ,GAMU) , 

PH2 

0450 

14Z( 12),GAMX), 

(2(13) »ETH) , 

(2114), FFA) , 

(2(15) .FFB) , 

PH2 

0460 

A( Z ( 16) , TMD2 ) , 

(2(17) ,IMX2) , 

(2(18), XMAX), 

(2(19) , TAM AX) , 

Ph* 

0470 

5(2(20), TYMAX)  , 

(2(21) ,AMOM ) v 

(2(22), AMXM), 

(2(23) ,ONN) • 

PH2 

C48C 

6(2124) ,OMINJ « 

(2(25)  .FEF)  , 

(2(26), OTNA) « 

(2(27 ) yCVIS)  , 

PH2 

C49C 

7(2(28) ,NPR) # 

(2(29) ,NPR1) , 

(2(30), NO, 

(2(31), NPC ) , 

PH2 

0500 

812(32) ,NRC) , 

(2(33), IMAX), 

12(34), JMAXA), 

(2(33) ,JMAX)  , 

PH* 

0510 

9(2(36), JMAXA), 

(2(37) , XMAX) • 

(2(38.*  -  KMXA)  , 

(2(39) ,NMAX ) 

PH* 

0520 

GtUUl VALENCE 

( 2( AO) , NO) , 

(2(41), KUT), 

(2(A2) , I XMAX) , 

PH2 

0330 

1(2(43), NOD)  , 

( 2 ( AA) ,NOPR)  , 

(2(A5),N1MAX) , 

(2(46) ,NJMAX) , 

PH2 

C54C 

2(2147 ), 111 « 

(2(46) ,12)  , 

( 2(A9  )  ,1 3) , 

(2(30) , 1  A) , 

PH2 

0550 

3(2(51), Nl), 

(2(52), N2) , 

(2(53), N3), 

(2 (5  A) ,NA) , 

PH2 

0360 

4 


o n n o n  o  o  ooonnoonooo 
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L 

i 

L 


L 

L 

L 

t 

t 


t 

U 

L 

l  c 

l 

l 

r* 

U 


442455), N5), 

5( 2( 59) ,N9)  , 
6(2463) , TRAD) * 
742467)  , RADER)  t 
BTzrrrrjREZFcri, 
94Z475),TQ20NE), 
OEQUI VALENCE 
142(82) «CABLN) « 
242486) ,  WSGD) « 
3(2490) ,  Sli« 

4424  94) , S5 }  , 
5(2(98), S9), 


(2456) ?N6) , 
(2(60), N10), 
42464)  »XNRG )  • 
(2(68) ,RA0£T ) « 
C  Z4  72) ,RSTQP) • 
4Z4  76) ,ECK) * 
42(79), X2), 
(2(83),VISC), 
(2487) , WSGX) . 

( 2  491) , S2) , 
42(95), S6I, 


42457)*N7),  42(58), N8), 

( Z4  6 1 )  «N1U •  (2(6i) ,NRM) , 

42(65), 5N),  42(66), OXN), 

42(69), RA0E8) ,  42470) , OTRAO) , 

(2473), SHELL),  42(74) .BBOUNO) 
(2(77), SBOUNO) ,  (2478), XI) 

42480) ,Y1) «  4Z481),Y2), 

(2484), T),  (2485), GMAX), 

(2488), GMA OR),  42(89) ,GMAXR) , 
42(92) , S3) ,  4Z(93) ,S4), 

42(96), S7),  (2497), S8), 


(2499) ,S10) 

EQUIVALENCE) 24 127) , SSI), 424 12 8) ,SS2) ,4  2  4129) , SS3) 
1424 131) *SS5) , 424 132) , SS6) ,424 133) , SS7) • 42( 134) ,SS 
2424 136) ,SS10) , (2( 137) ,SS11) 

EQUIVALENCE) 24 140) , SNR) , 4  24 141) , STR) , ( 2  4142 ) , SNT) 
14 24 143) ,ST7) ,42 4 144), SNB), (24 145), STB) 


•(2(130), SS4), 

8)  •  4  24  135)  •  SS9)  • 


©EQUIVALENCE  4XX(2),X(1)),  4UR, UL, FLEET) ,  (UR4 100) , YAMC) 

14 PR( 100) , SiGC) , (PR, PL, GAMC) ,4UR,TA8), 

24UR( 161 ,AHK) ,  4UR(31),PK),  4UR(46),QK),  4YY(2),Y(1)) 

EQUIVALENCE) AM, ASN ),4AM(53),AST),4XL,ASNB), 

1 4  XL4  53) ,ASTB) , ( YL,RSN) , 4 YL453) ,RST ) • 4 IW1,S IG33) , 

24IMK53)  ,OUDOT)  , 4 1 W2,DVDQT )  ,  4 1M24 53)  ,0A1X) 

EQUI VALENCE (L,KL) ,4M,KAR), 4MA,KA), (MB, KALI , 

14  MC,KBR) , (MO, KB! , ( ME ,KBL) 

EQUIVALENCE  (24  105)  ,SNI.)  ,  424 106), STL) 


******  NOTE  1  MATERIAL  ONLY  4X))  ************ 

INPUT  READS  OIL  DUMP  TAPE  OR 
MILL  CALL  SUBROUTINE  SET'UP  WHICH 

WILL  MAKE  A  OUMP  TAPE  FOR  CERTAIN  TYPES  OF  PROBLEM 
(SEE  SECTION  ON  SET’UP) 

ALSO  CALCULATES  DX  AND  DY  AND  EQUATION  OF  STATE  DATA 
CALL  INPUT 

CDT  ROUTINE  CALCULATES  OT 4 HYDRO  TIME  STEP) 

AND  PRESSURES,  ADVANCE  CYCLE  NO.  ETC. 

10  CALL  CDT 

IN  EDIT,  DETERMINE  WHETHER  TO  EXECUTE  A  LONG 
•  PRINT,  A  SHORT  PRINT,  A  TAPE  OUMP,  ETC.  AND 
CALCULATE  TOTAL  ENERGY  IN  SYSTEM4C0MPARE 
M, TH  ETH)  TOTAL  MASS,  INTEGRATE  TOTAL 
COMPONENTS  OF  MOMENTA. 

CALL  EDIT 

CALL  Si ITET ( 1 »KOOOFX } 

C  SENSE  LITE  1  SIGNIFIES  THIS 


PH2  0570 
PH2  0580 
PH2  0590 
PH2  0600 
PH2  0610 
PH2  0620 
PH2  0630 
PH2  0640 
PH2  0650 
PH2  0660 
PH2  0670 
PH2  0680 


PH2  0690 
PH2  0700 
PH2  0710 
PH2  0720' 


PH2  0730 
PH2  0740 
PH2  0760 
MAIN0020 
MAIN0030 
MAIN0050 


MA1N0060 


MAIN0070 


MAIN0080 

MAIN0090 
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C  IS  THE  LAST  CYCLE  OF  THIS  RUN  $$$$$$$$$$$$$$$ 

C  LITE  TURNED  ON  IN  THE  EDIT  ROUTINE  ****** 

GO  TG( 30,20) .KOOOFX 

I  C  PHI,  INTEGRATE  THE  MOMENTA  EQS.  INTEGRATE 
C  ENERGY  EQUATIQNIONLY  CHANGES  DUE  TO  WORK 
C  TERMS).  NO  MOVEMENT  OF  NASS  HERE 

20  CALL  PHI 

c  ******  PH3  calculates  the  change  in  the  velocity 
C  COMPONENTS  AND  INTERNAL  ENERGY  DUE  TO  THE 
C  STRESSES  ACTING  ON  THE  CELL  ... 

CALL  PH 3 

C  TRANSPORT  MASS  ACROSS  BOUNDARIES  I  SOLVE 
C  MASS  TRANSPORT  EQ.)  TRANSPORT  TERMS  IN 
C  THE  MOMENTA  AND  ENERGY  EQS.  LEFT  OUT  OF 
C  PHI,  HERE  APROXIMATED  BY  MASS  MOVEMENT.  CONSERVE 
C  MASS*  MOMENTA  AND  TOTAL  ENERGY. 

CALL  PH2 
C 

I  C 

SO  TO  10 
30  CALL  EXIT 
END 


£ 
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MAIN010 


MAINOIK 


MAIN012( 

MAIN0131 

MAIN014( 

MAIN015( 

MAIN016( 

MAIN017C 


%  I8FTC  CAROS  LIST »OfcCK,R£f 
SUBROUTINE  CAROS 

DIMENSION  TABLE! 11 .CARD! 7) , CABLE!  i i 
COMMON  TABLE 

C  A  2  IN  COLUMN  1,  ROUTINE  WILL  FIX  THE 
C  FLOATING  PT.  NO* 

C  A  1  IN  COLUMN  lf  MEANS  THIS  IS  LAST  CARO  TO 
C  READ  IN. 

EQUI VALENCE ( TABLE  !  1 ! , CABLE ( 1 ) ! 

WRITE  <6* iO) 

1  READ  <5,  ilHENO#LOC.  NUMWPC*  (CARO!  I  i , 1= 1 .NUMWPC) 
WRITE  !6* 12 ) I END.LOC . NUMWPC* ! CARO! I ) • Is I* NUMWPC) 
DO  4  1*1, NUMWPC 

J=LOC+I-l 
JFC IEND-2J2.5.2 
5  CABLE! J )s IF  I X( CARO(  I  i  ! 

GO  TO  4 

2  TABLE! J )=CARD(I ) 

4  CONTINUE 

IFti£NO-li 1»3» 1 

3  RETURN 

C  FORMATS 

iO  FORMAT (20H1  RFM  INPUT  CARDS///! 

U  FORMAT!llti5»IIf 0P7E9.4) 

12  FORMAT! 1H  I4,I7,I3,IP7E14.6) 

ENO 
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CARD0010 

CARD0020 

CAR00030 


CAR00050 
CAR00070 
CAR00080 
CAR00090 
CAROOIOO 
CAROOllO 
CAR00120 
CARD0130 
CAR00140 
CAR00150 
CARGO 160 
CAROOl 70 
CARGO!  80' 
CAR00190 

CAR00210- 

CARG0220 

CAR00230 


. $ 1BFTC  SfcTUP  UST»DECK,REF 
SUBROUTINE  SETUP 

C  WILL  UNLY  GENERATE  (11  MATERIAL. 

C  PACKAGES  MUST  3E  RECTANGLES. 

C  ASSUMPTION  Of  =  DX  AND  =  DY 
C  LOAD  PK ( 4) = 1. 

M=PK( 4) 

C  LOAD  PK( 5) =RIGHT  BOUNDARY  OF  PELLET ( 11. 
MA=PK<5) 

I  C  LOAD  PK16)=BQTT0MI JJH  OF  PELLET. 

MB=PK( 61 

i  C  LOAD  PK(7)=TQP(J)  OF  PELLET. 

MC=PK( 7 1 

|  C  LOAD  PK ( 8)  =  1  • 

MD=PK(8) 

I  C  LOAD  PK(9)=RIGHT( II  BOUNDARY  OF  TARGET. 
ME=PK<91 

I  C  LOAD  PKI 10 1= BOTTOM! J 1  *1  Of  TARGET. 
MZ=PK(10) 

I  C  LOAD  PK( 11 1=TQP( J10F  TARGET. 

N=PK ( 11) 

C  LOAD  INITIAL  DENSITY  INTO  ZUlli. 

W  RHO=Z (111) 

1  C  LOAD  INITIAL  PELLET  VELOCITY  INTO  2(112). 

VTEF=Z( 112) 

-  KMAX= IMAX* JMAX+1 

KMAXA=KMAX+1 
JMAXA= JMAX+1 
IMAXA=IMAX*1 

C  CLEAR  ALL  CELL  ARRAYS. 

DO  1  K=1«KMAX 
U( K) =0.0 
V(K)=0.0 
P(K)=0.0 
AMX(K 1=0.0 
A|X(K)=0.0 

1  CONTINUE 
DX( 1)=DX(  1 ) 

XI 1)=DX4 1) 

WS=X(1)**2 
If IGAM) 4(2«  4 

4  PIDY=1. 

TAUI 1 )  =  DX( 1 ) 

GO  TO  3 

2  PIDY=3. 1415927 
TAU(11=WS*PIDY 

C  CALCULATE  DX,X,TAU 

3  DO  10  I =2» IMAX 

I  X( Ii  =  X( I-ll+DXI  1) 

DX( I )=DX( 1) 


f 


I 
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SETUOOK 


SETU098v 

SETU0991 

SETUIOOC 

SETU101C 

SETU102C 

SETU103G 

SETU104C 

SETU105C 

SETU106C 

SETU107C 

SETU108C 

SETU109C 

SETUllOO 

SETU111C 

SETU1120 

SETU1130 

SETU114G 

SETU115G 

SETU116G 

SETU117G 

SETU1180 

SETU119G 

SETU120C 

SETU121C 

SETU122G 


SETU1240 


SETU1260 

SETU127G 


MSA*X41)**2 
IF(GAM) 5* 6,5 

5  TAU4IJ  =  DXU> 

GO  TO  10 

6  TAU(  I)*PI0Y*4WSA-WS) 

MS-MSA 

10  CONTINUE 

yiu^oyii) 

C  CALCULATE  OY  AND  Y. 

00  20  J*2,JMAX 
Y4J)*Y( J-1I+DY41) 

OY4J)*OY41J 
20  CONTINUE 
E  T  H*  0  -  0 
00  30  I=M,MA 
K*<  M8— 1  J*IMAX*-I*i 

C  CALCULATE  MASS*  AND  VELOCITY  OF  PELLET. 
00  40  J«M8»MC 
AMXi K)=AHG*DY ( J )*TAU( 1 J 
V(K)»VT £F 

C  CALCULATE  TOTAL  ENERGY  (ETH.l 

ETH=ETH*-AMX(K4*C  V(  Kl **21/2.0 
40  K*K+IMAX 
30  CONTINUE 

C  CALCULATE  MASS  OF  TARGET. 

00  50  I=MD»  ME 
K*(MZ-11*IMAX*H-A 
DO  60  J=MZ  *  N 
AMX4 K1»RH0*DY4 J ) *T AU <  I 1 
60  K-K+IMAX 
50  CONTINUE 
IMAXS1MAX 
JMAX*JMAX 
SHELL=2.0 
CYCLE=0.0 
DT=O.G 
NMAX=0 
Nl=2 
N2*3 
N3=0 
N4=127 
XMAX*X< IMAX) 

TXMAX=XMAX*2.0 
YMAX=Y4 JMAX) 

TYMAXSSYMAX*2.0 
REWIND  7 
WS=555.0 

C  WRITE  OUTPUT  FOR  OIL  ON  TAPE. 

WRITE  (  7)WStCYCLE«N3 
WRITE  (  7J4Z4II ,1=1* 150) 
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SETU1280 


SETU1300 

SETU1310 

SFTU1320 


WRITE  (  7HUUJ ,  V(I>  ,AMX(IJ  ,  AiX(I)  ,PU)  rl=l.KMAXAI 
WRITE  I  7JX(0)v(X(l)fTAU(IivI-lfINAX) 

WRITE  (  7) I Y( I i t I  =  Of  JMAX) 

WS=666.0 

WRITE  I  7} WS  t  WS « WS 
REWIND  7 
RETURN 
END 


SETU1330 

SETU1340 

SETU1350 

SETU1360 

SETU1370 

SETU1380 

SETU1390 


SETU1400 

SETUI410 

SETU1420 


SETU1430 

SETUI440f 

SETU1450 

SETUI460- 

SETUI470 

SETUI480 

SETU1490 

SETU1500 

SETU1510 

SETU1520 

SETUI530 

SETU1540 

SETU1350 

SETUI560 

SETU1570 

SETU1580 

SETU1590 

SETU1600 

SETUI610 

SETUI620 

SETU1630 

SETU1640 

SETU1650 


SETUI670 


o  o  o  o  n  n  no  no  non  noon  on  n  on  non 


TIBFTC  INPUT  LIST, DECK, REF 


SUBROUTINE  INPUT 

TURN  ON  SENSE  LITE  3. 

I NPU0010 
INPU0760 
INPU0900 

CALL  SLITE  (3) 

I NPU0980 
'NPU0990 

READ  HEADER  CARD  (COLUMNS  2-72). 

READ  ( 5, 8004J IMS 

INPUXOOO 

WRITE  (6, 800411**$ 

CALL  DATA. 

INPU1010 

6 

CALL  CARDS 

IF  PK13J  *  OR  GREATER  THAN  ZERO,  CALL  ROUTINE 

SET-UP,  OTHERWISE,  BINARY  OIL  TAPE  HAS  BEEN  MADE. 

INPU1020 

READ  IN  DATA  FROM  OIL  DUMP  TAPE,  OR 

GENERATE  A  DUMP  TAPE  FOR  OIL,  AND 

CALCULATE  OX  AND  OY  FROM  THE  X  AND 

Y  VALUES  FROM  TAPE. 

IF C  PM  31 18887, 8888, 8888 

INPU1030 

6888 

CALL  CARDS 

INPU1040 

CALL  SETUP 

INPU1050 

8887 

CONTINUE 

INPU1060 

INPU1070" 

READ  TAPE 

GO  READ  BINARY  TAPE. 

INPU1080 

GO  TO  1000 

INPU1090- 

INPU1100 

READ  IN  REMAINING  INPUT  CARDS 

INPU1110 

10 

CONTINUE 

INPU1120 

CALL  CARDS 

INPU1130 

GO  TO  2000 

SET  THE  PRESSURES  TO  ZERO. 

INPUL140 

INPU1150 

40 

DO  45  K~1«KMAXA 

INPU1160 

45 

P(K)=0.0 

INTEGRATE  BACKWARDS  ON  CYCLE,  TIME  ANO  NO.  OF 

CYCLES  BETWEEN  ENERGY  CHECK,  SINCE  THESE 

ARE  ADVANCED  IN  CDT • 

NOTE,  RSTOP  =  ENERGY  FACTOR  AND 

SHELL  =  FACTOR  ON  SPEED  OF  SOUND  CALC. 

CONVERT  FROM  JERKS  TO  ERGS. 

Z(117)=Z(117)*RST0P 

ZI1191=Z(119)*RST0P 

Z(1221=Z(122)*RSTGP 

Zt 126J=Z( 126)  *RSTQP 

Z( 110  )=Z( 110)*RST0P 

INPU1170 

RST0P=1.0 

INPU1178 

T =T-DTNA 

INPU1180 

NC=NC-1 

IF (Z( 149)13000,3001 ,3000 

INPU1 190* 

ooooo  o  o  o  o  r>  r>  no 
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3000  Z (148 1 =S0RT ( Zl 119) /Z (1151) 

Z ( 149 ) =Z ( 149) *SHELL 
SHELL=1.0 

3001  C YCLE-NC 
NPC=NPC-1 
UVMLX=0.G 

CALCULATE  THE  DX«S,  SINCE  THESE  ARE 
TAPE. 

DO  50  1  =  1, IMAX 
50  DX(I)=XiIl-X(I-l) 

CALCULATE  THE  DY»S,  SINCE  THESE  ARE 
TAPE. 

DO  55  J=l»  JMAX 
55  DYI J)=Y(J)-Yl J— 1 ) 

J=MZ— 8 

PRINT  Z  BLOCK. 

62  DO  80  1=1, J,8 
K= I  +7 

DO  65  J=  I » K 
IF(Z( J))70,65,70 
65  CONTINUE 
GO  TO  80 
70  K=I+7 

WRITE  (6,811111, (2(L) ,L= I , K) 

80  CONTINUE 


INPU12QL 

INPU12K 
I NPU122C 

NOT  ON 

1NPU123C 
INPU124C 
ON 

INPU1251 
INPU1261 

y 


I NPU129C 
INPU130C 
INPU131C 
I NPU132C 
INPU1331 
INPU134C 
INPU135C 
INPU136C 


ASSUMPTION  THAT  ALL  DX  AND  OY  ARE  = 
NOTE,  DVK=K(0J/IRHO(0)*DX  SQ.l 
WS=DX ( 1 ) *DX( 1 ) 

DVK=DDXN/(ZI 111J*WS) 


GO  TO  10000 

READ  BINARY  TAPE. 

INPUI37 

INPU138 

INPU139 

INPU140 

1000 

MZ=150 

INPU141 

I  WS=0 

INPU142 

1003 

REWIND  7 

1004 

READ! 7) PRI 11 ,PR(2) ,N3 

NR=N3+5 

INPU145 

1006 

IFIPRt 11-555.01 1010,1016,1010 

INPU146 

1010 

IWS=IWS+1 

INPU147 

1011 

I F ( MUD ( IWS, 3119902, 9902, 1003 

INPU148 

1016 

IF( PR (2)) 10 10,1018, 1018 

CHECK  HERE  FOR  THE  CORRECT  CYCLE  NUMBER. 

INPU149 

1018 

IF (PK(21-PR( 2)) 1023,1023,1020 

SKIP  OVER,  LOOK  AT  NEXT  CYCLE. 

INPU150 

1020 

DO  1P^2  L=2 *NR 

INPU15 I 

1022 

R£AO( 71 
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GO  TO  1004 

INPU153Q 

1023 

READ! 711Z1 1 1  »i=l»MZ) 

C 

CHECK  FOR  THE  CORRECT  PROBLEM  NO. 

1 F4 ABS(PROB-PK( 1) l-.Oll 1024, 1024,9901 

INPU1550 

1024 

READ!  7)  1U<  I1..VU1  ,AMX(I)  ,AIX(I),P1I  ),i=i,KMAXA) 

REAO( 7 1 X ( 01 , (X( 1 1 ,TAU( I 1 , 1=1 , IMAX) 

REAOI 71  ( Y(  1 1 ,  1=0, JMAX) 

1025 

CONTINUE 

INPU1650 

1034 

READ(7IPR( 1), PR (21, PR (3) 

1036 

1F1PRUJ-555.019904, 1040, 1038 

INPU1680 

1038 

IF (PR( 21-666. 0)9905, 1040, 9905 

1NPU1690 

1040 

GO  TO  10 

1NPU1700 

c **** 

END  OF  READ  TAPE  *************************«***********************INPU1710 

C 

I NPU1720 

C 

INPU1730 

c 

CALCULATE  MAX.  GAMMA  AND  GAMMA/ ( GAMMA- 1.1. 

c 

1 NPU1740 

2000 

I F ( WSGX 19906,2010,2005 

INPU1750 

2005 

GAMX=1.0/< WSGX-1.01 

1NPU1760 

2010 

WSGX=<GAMX«-1. 0)  /GAMX 

INPU1770 

GMAXR=GAMX*MSGX 

1NPU1780 

2012 

I F( WSGDT 9907, 2020, 2015 

INPU1790 

2015 

GAMD=1. 0/( HSGD— loOJ 

INPU1800 

2020 

WSGD=i GAMD+1 .0 ) /GAMD 

INPU1810 

GMADR=GAMD*NSGD 

I NPU1820 

GMAX=WSGD 

INPU1830 

1F( WSGO-WSGX 12025,2030, 2030 

1 NPU1840 

2025 

GMAX=WSGX 

I NPU1850 

2030 

GO  TO  40 

INPU1860 

c«*** 

END  OF  R  E  S  ♦♦***************************************************(NPU1870 

C 

1NPU1880 

C 

INPU1890 

C 

ERROR 

INPU1900 

9901 

NK=1023 

INPU1910 

GO  TO  9999 

1 NPU1920 

9902 

NK=1011 

INPU1930 

GO  TO  9999 

INPU1940 

9904 

NK=1036 

INPU1950 

GO  TO  9999 

INPU1960 

9905 

NK=1038 

INPU1970 

GO  TO  9999 

INPU1980 

9906 

NK=2000 

INPU1990 

GO  TO  9999 

INPU2000 

9907 

NK=2012 

INPU2010 

9999 

NR=1 

I NPU2020 

CALL  OUMP 

INPU2030 

C 

1NPU2040 

10000 

RETURN 

INPU2050 

C 

1NPU2060 

C 

FORMATS 

INPU2070 

I 


o n o o n  o  o  o  o  o  ooooooooo 
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pS  $  IBFTC  COT  LIST, DECK, REF 
i  SUBROUTINE  COT 


CHECK  COURANT  CONDI TiON  AND  PARTICLE 
VELOCITY. 

RECORD  I  ANO  J  OF  ZONE  WHERE  OT  IS  BEING 
CONTROLLED. 

3000  VEL=0.0 
3005  DO  3050  1=1,11 
3010  K=I*1 
3015  DO  3050  J=1,I2 
1  =  1 
J=J 

3020  IF  I AMXIK) i 9901,3050,3025 

CALCULATE  PRESSURES  FROM  EQUATION  OF  STATE (ESI. 
3025  CALL  ES 

3030  IF  I ASS( P(K) )— 1.0E—20 13035,3035,3040 
3035  P( K) =0. 0 

3040  IFIWSGX-VED3050, 3050,3045 
3045  VEL*WSGX 
3050  K=K*IMAX 
3055  KDT=1 

UVMAX=-i.O 
3070  DO  3255  1=1,11 
3075  K*I*1 
3095  DO  3255  J=1,I2 
3100  KP=K*IMAX 

IF (AMXIK) 19901 ,3255,4 

IF  RH01K)  IS  LESS  THAN  2(138),  CELL  K 

MILL  BE  BYPASSED  FOR  STABILITY  CHECK. 

4  IF (ANX( KJ / I TAUII1*DY!J) 1-2(138) )3255, 3255, 3115 
3115  SIG=DX(I) 

3120  IF  lDYUJ-SIGm25,3130,3130 
3125  S1G=CY(J) 

C=SPEED  OF  SOUND  FOR  POLYTROPIC  GAS  AS 
THE  SQ.  ROOT  OF  I GAMMA*P/RHO) . 

HERE  CALCULATE  THE  SPEEC  OF  SOUND  FOR 
THE  EQUATION  OF  STATE 
AS  THE  SQ.  ROOT  OF  DP/DRHO. 

3130  IF(Z( 148) 14000,4000,4001 

4000  WS=SQRT(GMAX*TAU(I)*DY( J)*ABS( PIK) l/( AMX(K) ) ) 

GO  TO  3205 

4001  WSA=A8S ( PI KJ ) 

WS=Z ( 148)+Z( 149) *<  WSA**ZI 150) ) 


CDT  0010 
CDT  0020 
CDT  0990 
CDT  1000 
CDT  1010 
CDT  1020 


CDT  1030 
CDT  1040 
CDT  1050 
CDT  1060 
CDT  1070 
CDT  1080 
CDT  1090 
CDT  1100 

CDT  1110 
CDT  1120* 
CDT  1130 
CDT  1140 
CDT  1150. 
CDT  1160 
CDT  1170 
CDT  1180 
CDT  1190 
CDT  1200 
CDT  1210 
CDT  1220 
CDT  1230 
CDT  1240 


CDT  1250 
CDT  1260 
CDT  1270 
CDT  1280 


CDT  1290 
CDT  1300 
CDT  1310. 

CDT  1330 
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3205  WS=WS/SIG 

3210  IF (UVMAX-WSI3215, 3220*3220 
3215  N10=I 
Ni  1=J 
UVMAX=WS 
3220  CONTINUE 

C  EULER  I AN  CHECK  FOR  RADIAL  PARTICLE  VELOCITY. 

1  IF  ( GAM)  2)3)2 

3  WS“ABSIU1K))/TAU(I)*X(I)/.5*PI0Y 
GO  TO  3225 

C  FOR  CARTESIAN  CODE 

2  WS=ABS(UIK) )/DX(  I) 

3225  IF ( UVMAX-WS 13230*3235  *3235 
3230  UVMAX=WS 
N 10=  I 
N 1 1= J 

3235  WS-ABS (VlKii/DYlJI 
<3240  IF( UVMAX— WSJ 3245 *325 0*32 50 
3245  N10=I 
Nl  1- J 
UVMAX=WS 
3250  CONTINUE 
3255  K=K+ I MAX 

IF (UVMAXl 9912* 9912 *3260 
C  FOR  OPTIONS  ON  CABLN,  CHECK 

C  SECTION  3.4  IN  GAMD-5580. 

3260  IF(CABLN)90, 91,3300 

90  0T“. 5/VEL/UVMAX*Z 1 139 J 
GO  TO  3295 

91  WS=UVMAX*DT 
WSA=0.5/VEL 

3265  IF (FF A- WSAJ  3276*32  76*3270 
3270  FFA=WSA 

3276  IF(WS-FFAJ3285*3300,3280 
3280  0T-DT/RS*FFB/0.9 
GO  TO  3295 

3285  I F< WS— FFB 13290*3290 » 3300 
3290  DT=DT*FFA/WS*0.9 
3295  KDT=0 

C  INTEGRATE  THE  TIME  AND  CYCLE  COUNTER. 

3300  T= J+DTNA 

85  IF (DTRA 0)9911*80 *81 

80  NR=NRM 
84  WS=NR 

TRAO=DT/WS 
GO  TO  82 

81  I WS=DT/DTRAD 
NR=IMS+1 

83  IFINR—NRM) 84,84,80 

82  NC=NC+1 


CDT  1350 
COT  1360 
CDT  1370 
CDT  1380 
'  l  CDT  1390 

0 


COT  1420 
COT  1430 

CDT  1440 
COT  1450 
CDT  1460 
CDT  1470 
COT  1480 
CDT  1490 
CDT  1500 
CDT  1510 
COT  1520 
CDT  1530 
CDT  1540 
CDT  1550 


CDT  1560 
CDT  1570 
CDT  1580 
CDT  1590 
CDT  1600 
CDT  1610 
CDT  1620 
CDT  1630 
CDT  1640 
CDT  1650 
CDT  1660 
CDT  1670 
CDT  1680 

CDT  1690 
CDT  1700 
CDT  1710 
CDT  1720 
CDT  1730 
CDT  1740 
CDT  1750 
CDT  1760 
CDT  1770 
CDT  1780 
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CYCLE=NC 

NPC-NPC+1 

3305  If 11)9909,3320,3310 
3310  If (KDT)9910, 3315, 3320 
3315  WRITE  16, 8000)1, OTNA.DT 
3320  DTNA=DT 

GO  TO  3325 

C  NEGATIVE  NASS 

9901  NK=3020 

GO  TO  9999 

9909  NK=3305 

GO  TO  9999 

9910  NK=3310 

GO  TO' 9999 

C  THE  OT  WILE  BE  0.  OR  NEGATIVE  .STOP 
9912  NK=1 

GO  TO  9999 

9911  NK=85 
9999  NR=2 

CALL  DUMP 
3325  RETURN 

SOOOOf ORMAT  I 17H0CHANGE  DT  ...  T=1PE9.3,UH 
II N^I)=1PE9«3) 

END 


CDT  1790* 
CDT  1800 
CDT  1810 
,  \  CDT  1820 

CDT  1830 
CDT  1840 
CDT  1850 
CDT  1860 
CDT  1870 
CDT  1880 
CDT  1890 
CDT  1900 
CDT  1910 
CDT  1920 


CDT  1930 

CDT  k  «•* 
CL/  < 

CDT  lv6G 

OT (N)= 1PE9.3, 13H  DTCDT  1970 

CDT  1980 
CDT  199C 


onnnn on  n  no  on  oooooooooooo 
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$1 BFTC  PHI  LIST, DECK, REF 

SUBROUTINE  PHI 

VELOCITIES,  ENERGIES,  PRESSURES  ARE  AT  THE 
CENTER  OF  THE  CELL. 

(2i  PASSES  THRU  PHI  ARE  REQUIRED.  NO 
MASS  IS  MOVED  IN  PHI. 

******  NOTE  1  MATERIAL  ONLY  (XI J  ************ 


NRT =0 
NRC=0 
UU=1®  E+ 15 
UT=0.0 

YOU  WILL  GET  BACK  HERE  IF  AIX  WAS  LESS 
THAN  0.  AND  PROVIDED  SN=0. 

8000  VEL=1 .0 

INITIALIZE  MID-POINTS  OF  FIRST  AND  SECOND 
CELL  IN  R  DIRECTION. 

1F(GAM}9000, 3301, 9000 
9000  RC-1. 

RR-RC 

GO  TO  3304 

3301  RC=DX(l)/2.0 
RR=(XUJ*X12)  J/2.0 

3304  K=2 

AXIS  OF  SYMMETRY  BOUNDARY  CONDITIONS. 

DO  3302  J= 1 »  JMAX 
PL< JI=PIKJ 
UL  ( J ) =0*0 

3302  K=K+ I MAX 

FIRST  PASS  THRU,  CALCULATE  U  AND  V  AT 
CYCLE  NU,  AND  THE  WORK  TERMS  USING  U  AND  V 
FROM  CYCLE  N„ 

SECOND  PASS  THRU,  CALCULATE  ONLY  THE 
CONTRIBUTION  TO  THE  CHANGE  IN  INTERNAL  ENERGY 
FROM  WORK  TERMS  EVALUATED  FROM  U  AND  V 
AT  CYCLE  N+ 1. 

DO  3360  1=1,11 
K=  I  + 1 

IF (C VIS J 7002, 7003, 7003 
C  BOTTOM  BOUNDARY  IS  TRANSMI TT I VE . 

7002  VBLQ=  VI K ) 

PBL0=0. 0 
GO  TO  7004 

C  BOTTOM  BOUNDARY  IS  REFLECTIVE. 


PHI  001 
PHI  090 

( 

( 

PHI  09b 
PHI  099 
PHI  100 
PHI  101  ( 

PHI  10?  ‘ 

PHI  103- 
PH1  104 
PHI  105 
PHI  106 
PHI  107  1 

PHI  108  { 


PHI  109 

< 

I 


PHI  110 
PHI  111 
PHI  112 

PHI  115 
PHI  114 
PHI  115 
PHI  life 


I 

i 

PHI  117 
PHI  118 
PHI  115  ( 

PHI  120 
PHI  121 ! 

PHI  122  ( 

I 

\ 


no  o  no  n  on  on  on  non  on  non 
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7003  VBLO=0.0 
PBLO=P(K) 

7004  TAUDTS=TAU( I )  *DT 

11=  MAX*  I I )  OF  DISTURBANCE  IN  R  DIRECTION. 

12=  MAXI J )  OF  DISTURBANCE  IN  2  DIRECTION* 

00  LOOP  IN  J  DIRECTION 
DO  3348  J=1 « 12 
PIDTS=1.0/(PIDY*DT*DYI JJ J 
IFIGAMi 9002 , 9004, 9002 
9002  P IDT S=2 **PIDT S 

K=  INDEX  OF  CELL  IN  QUESTION. 

N=  INDEX  OF  CELL  ABOVE. 

9004  N=K*IMAX 

3305  IFIAMXIK) >9902,3340,3306 

3306  IF( IMAX-IJ 9903*331 1 *3310 

3310  IF(AMX1K+1)19904, 3312, 3314 

ME  ARE  AT  THE  RIGHT  BOUNDARY,  SET  PRESSURE 
GRADIENT  TO  0.  IN  R  DIRECTION,  MODIFY  ETH. 

FOR  RIGHT  BOUNDARY  BEING  TRANSMITTI VE. 

3311  PRR=PLIJJ 

3307  £TH=EFH-PRR*UIK)/PIOTS*RC 
GO  TO  3313 

RIGHT  BOUNDARY  CONDITION  FOR  THE  MOMENTUM  EQ. 

ADJACENT  TO  EMPTY  CELL. 

3312  PRR=0. 0 

3313  URR=RC*U(KJ 
GO  TO  3316 

CALCULATE  PRESSURE  AT^lNTERFACE 1 II  AND 
IRUI  FOR  WORK  TERM. 

3314  PRR=(PI  Kl+P ('{■*■  U  1/2.0 

3315  URR=(U(K)*RC*U(K+lJ*RR)/2.0 

3316  IFIJMAX-J19905, 3318,3320 

SET  PRESSURE  GRADIENT  TO  0.  THIS  IS  FOR  TOP 
BOUNDARY  BEING  TRANSMITTI VE. 

33 IB  PA80VE=PBL0 

MODIFY  ETH  FOR  TOP  BOUNDARY  CONDITION. 

3319  ETH=ETH-PABOVE*V(KJ/2.0*TAUOTS 
GO  TO  3323 

3320  IF! AMX(NJ 19906,3322, 3324 

TOP  BOUNDARY  CONDITION  I  EMPTY  CELL  ABOVE. I 

TOP  BOUNDARY  CONDITION  FOR  VELOCITY  I EMPTY  CELL  ABOVEI. 

3322  PABGVE=0.0 

3323  VABOVE=VI K ) 

GO  TO  3328 

CALCULATE  PRESSURE  AT  INTERFACE (Ji 

3324  PAB0VE=(Pm+P(N)!/2.0 
IFICVISi700i, 3325, 3325 

7001  I FI 1—J 1 3325*  7000*  9905 

BOTTOM  BOUNDARY  IS  TRANSMIT! IVE ,  SET  PRESSURE 
GRADIENT  TO  0. 


PHI 

1230 

PHI 

1240' 

PHI 

1250 

PHI 

1260 

PHI 

1270 

PHI 

1290 

PHI 

1300 

PHI 

1310 

PHI 

1320 

PHI 

1330 

PHI 

1340 

PHI 

1350 

PHI 

1360. 

PHI 

1370 

PHI 

1380 

PHI 

1390 

PHI 

1400 

PHI 

1410 

PHI 

1420 

PHI 

1430 

PHI 

1440 

PHI 

1450 

PHI 

1460 

PHI 

1470 

PHI 

1480 

PHI 

1490 

PHI 

1500 

rnn  o  o  n  n  o  o  no  no  no 
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.C  AND  MODIFY  ETH. 

7000  PBL0=PA80Vt  PHI  It 

ETH=ETH+PBLG*V(K)/2.0*TAUDTS  PHI  It 

C  VELOCITY  AT  I  NT  ERF ACE ( J ) 

3325  VABGVE-(VIK)+V(Nl)/2.0  PHI  If 

3328  IFIVEL19907, 3404, 3400  PHI  if 

C  COMPUTE  OELTA  U  AND  DELTA  V- 

3400  V(Kl=V(Kl+(PBLO-PABOVE)*TAUDTS/IAMX{Kl)  PHI  It 

****  NOTE,  EPSILON  IS  FOR  C.G.S.  UNITS  ************* 

***  FOR  OIL  UNITS  SET  IT  TO  l.E-8  ******* 

IFIABSi V(K)1~1. 13401,3401,3402 

3401  V*K1=0.0  PHI  1 

3402  UIK)=U(K)+(PL(J  J— PRR ) / ( AMX ( K J  1 *RC/P InTS*2. 0  PHI  If 

IFlABSlU(K) J-l. 13403, 3403,3404 

3403  U(K)=0.0  '  PHI  It 

CHECK  FOR  ADVANCING  COUNTERS  OF  THE  ACTIVE 
GRID  IN  THE  R  DIRECTION.  , 

3404  IF( 1-1116016,6005,6005  PHI  It 

6005  IFIUCK) 16605,6606,6605  PHI  it 

6605  NRC=1  PHI  If 

6606  I  FI  VC  Kl 16607 ,6004,6607  PHI  1( 

6607  NRC=1  .  PHI  1/ 

6004  IFCAIX(K)16015,6016,6015  PHI  1< 

6015  NRC=1  PHI  It 

6016  CONTINUE  ,  PHI  It 

HERE  CALCULATE  CHANGE  IN  INTERNAL  ENERGY 
DUE  TO  WORK  TERMS  ONLY. 

WS=  I VBLO-VABOVE ) *TAUDTS/2.0*P(K1  PHI  it 

RHQ=WS+ ( UL1 J )— URR) /P I DTS*P ( K J  PHI  1 

CONVERT  TO  SPECIFIC  INTERNAL  ENERGY. 

3332  WSX=A1X1K1 +RHQ/AMX(K)  PHI  1 

GO  TO  1000  PHI  1 

CHECK  FOR  NEGATIVE  INTERNAL  ENERGIES. 

1000  IFIWSXU011, 1001, 1001  PHI  1 

1001  A IX ( K 1 = WSX  PHI  1 

GO  TO  3342  PHI  1 

1011  UT= 1.0  PHI  1 

COMPUTE  NEW  DTI  STORE  IN  UU1  ASSUMING 
THAT  OI/DT  WILL  BE  THE  SAME  FOR  A  SMALLER 
TIME  STEP,  THE  NEW  DT  IS  CHOSEN  SUCH 
THAT  A I XI  AT  N+l 1  =  2/3  OF  AIXCN1. 

WSA=2.0*AI X  IK1/3.0*DT/(AIX(K1-WSX1  PHI  1 

1013  IF(WSA-UU) 1014,1001,1001  PHI  1 

1014  UU=W$A  PHI  1 

GO  TO  1001  PHI  1 

CELL  (KJ  IS  EMPTY,  SET  INTERFACE  QUANTITIES, 

ASSUMING  CELL  TO  THE  RIGHT  AND  TOP  ARE 
NOT  VOID. 

.  3340  PRR=0«G  PHI  1 

URR=UC K+l 1 *RR  PHI  it 
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PA8GV£=0.0 

VABGVE-VIN) 

C  SET  RIGHT  QUANTITIES  TO  THE  LEFT  IFOR  NEXT 
C  COLUMN  SHEEP  I  AND  SET  A8GVE  QUANTITIES  TO 
C  BELOW  FOR  NEXT  CELL  ABOVE. 

3342  V8L0=VA80VE 
PLI J»=PRR 
ULI J)=URR 
K=N 

3348  P8L0=PAB0Vr 
LL=K—l Max 

C  CHECK  FOR  ADVANCING  COUNTERS  OF  THE  ACTIVE 
C  GRID  IN  Z  DIRECTION. 

I  FI U ILL  J ) 6000 ,6001*  6000 

6000  NR  S= I 

6001  1F1 VILLI  16002*6003*6002 

6002  NRI=i 

6003  IFIAIXILLJ J6017, 6018, 6017 

6017  NRI=1 

6018  CONTINUE 
3355  RC=RR 

I FI GAM! 3360, 9007, 3360 
9007  RR=( XI I  +  IJ^XI 1+21  i/2.0 

3360  CONTINUE 

3361  IFIVELJ9911, 10000,3363 
3363  VEL-0.0 

GO  TO  3301 

C  ERROR 

9902  NK~^a05 

GO  TO  9999 

9903  NK=3306 

GO  TO  9999 

9904  NK=33I0 

GO  TG  9999 

9905  NK=3316 

GO  IC  9999 

9906  NK=332Q 

GO  TO  9999 

9907  NK-3328 

GO  TO  9999 
9911  NK=3361 
9999  NR=3 

CALL  DUMP 

C  IF  SNlftCT-0. J  ANY  NEGATIVE  ENERGIES  HILL 

C  REMAIN.  iF=0  »  CODE  WILL  TRY  ANOTHER  PASS 

C  WITH  A  SMALLER  DT. 

10000  IFISN17030, 7031,7030 
7031  IFIUT 17020*  7030*7010 
C  NEGATIVE  ENERGIES  HAVE  OCCURED,  INTEGRATE 

C  BACK  TO  CWCLE  N  WITH  l-DTU 


PHI  1830 
PHI  1840* 


PHI  1850 
PHI  1860 
PHI  1870 
PHI  1880 
PHI  1890 
PHI  1900 


PHI  1910 
PHI  1920 
PHI  1930 
PHI  1940 
PHI  1950 
PHI  1960 
PHI  1970 
PHI  1980 


PHI  2000 
PHI  2010 
PHI  2020. 
PHI  2030 
PHI  2040 
PHI  2050 
PHI  2060 
PHI  2070 
PHI  2080 
PHI  2090 
°H1  2100 
PHI  2110 
PHI  2120 
PHI  2130 
PHI  214C 
PHI  2150 
PHI  2160 
PHI  2170 
PHI  2180 
PHI  2190  - 


PF’I  2200 
PHI  2210 
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.  7010  UT=— 1.0 
DT=— DT 

YOU  NOW  HAVE  INTEGRATED  BACK  TO  CYCLE  N.  NOW 
INTEGRATE  TO  CYCLE  N*1  WITH  NEW  DT I STORED  IN  UU). 
GO  TO  8000 
7020  UT=0«  0 
DT=UU 

NR=DT/TRAD+ 1.0 
WS=NR 

TRAD=DT/wS 
DTNA=DT 
GO  TO  8000 

C°  INCREASE  ACTIVE  GRID  COUNTERS  IF  NEEDED. 

7030  11=1 1+NRC 
I2=I2+NRT 

IF (11- 1 MAX 16100,6100, 6200 
620C  1 1=1  MAX 

6100  IFI I2-JMAX16201, 6201, 6202 
6202  1 2  =  JMAX 
6201  RETURN 
END 


PHI  2220 
PHI  2230 


PHI  2240 
PHI  2250 
PHI  2260 
PHI  2270 
PHI  2280 
PHI  2290 
PHI  2300 
PHI  2310 

i 

PHI  2320 
PHI  2330 
PHI  2340 
t  PHI  2350 
PHI  2360 
PHI  2370 
PHI  2380 
PHI  2390 


uuuoeoouugouuuuuoooouuaoouuo 
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$IBFTC  PH 2  LIST  ,  DECK, REF 

SUBROUTINE  PH2  PH2  0010' 

NOTE  MIN.  DENSITY  FOR  REZONE  IS  A  INPUT  NO.  IVTII 
Z(U0)  =  CRITICAL  ENERGY  I  BETWEEN  GAS  ANO  CONDENSED  STATE! 

ZI1U)  =  INITIAL  DENSITY 
ZU12)  =  INITIAL  VELOCIJ/  OF  PELLET 
TOZONE  =  MINIMUM  DENSITY  FOR  MASS  FLOW 

PH2  0900 
PH2  0980 

AMPY=MASS  ACROSS  TOP  BOUNDARY  OF  CELL 
AMUT-RADIAL  MOMENTA  OF  THIS  MASS 
AMVT=AXIAL  MOMENTA  OF  THIS  MASS 
DELET=TOTAL  SPECIFIC  ENERGY  OF  THIS  MASS 

AMMP=MASS  ACROSS  RIGHT  BOUNDARY  OF  CELL 
AMUR=RADI AL  MOMENTA  OF  THIS  MASS 
AMVR=AXIAL  MOMENTA  OF  THIS  MASS 
OELER=TOTAL  SPECIFIC  ENERGY  QF  THIS  MASS 

AMMY=MASS  ACROSS  BOTTOM  BOUNOARY  OF  CELL 
AMMU=RADIAL  MOMENTA  OF  THIS  MASS 
AMMV=AXLAL  MOMENTA  OF  THIS  MASS 
D£LEB=TOTAL  SPECIFIC  ENERGY  QF  THIS  MASS 


GAMC=MASS  ACROSS  LEFT  BOUNDARY  OF  CELL 
FLEFT=RADIAL  MOMENTA  OF  THIS  MASS 
YAMC=AXIAL  MOMENTA  OF  THIS  MASS 
SIGC=TOTAL  SPECIFIC  ENERGY  OF  THIS  MASS 


NRT=0 
NRC=0 
REZ=0 .0 

CALL  SLITE  40) 

PIDTS=1.Q/ ( PIOY*DT I 

101  DO  103  J=1,JMAX 

102  GAMCI J)=0.0 
FLEFTI J 1=0.0 
YAMCI J!=0.0 
SIGCI Ji=0.0 

103  CONTINUE 

104  DO  547  1=1*11 
J=I 

105  K=I*1 

80  IFIAMX(KJ!9900*97,81 

81  IFIV(X) )82, 97,97 
97  AMMV=0.0 

GO  TO  98 

82  AMMY=AMXIK) *V(Ki*DT/DY<  JI 

83  IFI AMMY+AMXIK) J  84,85,85 


PH2  0990 
PH2  1010 
PH2  10 tO 
PH2  1030 
PH2  1040 
PH2  1050 
PH2  1060 
PH2  1070 
PH2  1080 
PH2  1090 
PH2  1100 
PH2  1110 
?H2  1120 
PH2  1130 
PH2  1140 
PH2  1150 
PH2  1160 
PH2  1170 
PH2  1180 
PH2  1190 
PH2  1200. 
PH2  1210 
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» 

64 

AMMY=— AMX ( K 1 

PH2 

1223 

85 

IF(CV1S)106,99,99 

PH2 

1230 

c 

BOTTOM  BOUNDARY  IS  TRANSMI TT I VE  ,  MATERIAL 

IS  MOVING 

c 

OUT,  RfcMQVE  ITS  ENERGY  FROM  ETH. 

106 

AMMU=AMMY*U(K) 

PH2 

124C 

AMMV=AMMY*VIK) 

PH2 

1250 

DElEB^AIXIK  )+(U(Kl**2+V(KI**2)/2.0 

PH2 

1260 

WS={UtKJ**2+VlK)**2)/2.0 

PH2 

1270 

£TH=ETH*AMMY#IAIXI  KJ+WS ) 

PH2 

1280 

GO  TO  107 

PH2 

1290 

c 

BOTTOM  BOUNDARY  IS  REFLECTIVE,  NET  MOMENTA 

CHANGE 

c 

IN  Z  DIRECTION  IS  2  MV. 

99 

AMMV=2.0*AMMY*VIK) 

PH2 

1300 

98 

AMMY=0.0 

PH2 

1310 

AMMU=0. 0 

PH2 

1320 

DELEB=0.0 

PH2 

1330 

c 

BEGIN  DO  LOOP  IN  J i Z )  DIRECTION. 

107 

00  546  J=1 , 12 

PH2 

1340 

108 

L=K+ IMAX 

PH2 

1350 

1  =  1 

PH2 

1360 

J=J 

PH2 

1370 

AREA=0.0 

PH2 

1380 

VEL=0 .0 

PH2 

1390 

FS=0.0 

PH2 

1400 

210 

IF( JMAX— J)211,211,207 

PH2 

1410 

4 

211 

VEL=1 .0 

PH2 

1420 

GO  TO  208 

PH2 

1430 

207 

IFIAMX(L) 1215,215,214 

PH2 

1440 

214 

IF(AMXIK) 1216,216,209 

PH2 

145C 

216 

VABOVE=V(U 

PH2 

1460 

GO  TO  212 

PH2 

1470 

215 

IFIAMX(K) 1205,205,208 

PH2 

1480 

205 

VABQVE=0.0 

PH2 

1490 

GO  TO  212 

PH2 

1500 

208 

VABQVE=V<K1 

PH2 

1510 

GO  TO  212 

PH2 

1520 

209 

VABQVE=(ViK!*V(LH/2.0 

PH2 

1530 

212 

CONTINUE 

PH2 

1540 

1  =  1 

PH2 

1550 

J  =  J 

PH2 

1560 

FS=0.0 

PH2 

1570 

404 

IFUMAX-1 1412,412,405 

'"'H2 

1580 

405 

IF (AMX1K+11 1411,411,409 

PH2 

1590 

409 

IF( AMXIK) 1410,410,407 

PH2 

1600 

410 

URR=U ( K+ll 

PH2 

1610 

GO  TO  408 

PH2 

1620 

411 

IF(AMX(K)1403,403,406 

PH2 

1630 

403 

URR=0.0 

PH2 

164o 

i 

f 

GO  TO  408 

PH2 

1650 

c 

WE  ARE  AT  THE  RIGHT  BOUNDARY  OF  THE  GRID, 

THE 

o  o 
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C  BOUNDARY  CONDITION  HERE  IS  TRANSMI III VE. 


412 

FS=1.0 

PH2 

1660' 

406 

URR*U(K) 

PH2 

1670 

GO  TO  403 

PH2 

1680 

407 

URR=(U(KJ+U(K+in/2.0 

PH2 

1690 

40B 

CONTINUE 

PH2 

1700 

109 

IF(AREAi990i,30i,547 

PH2 

1710 

301 

IF( VAB0VE1 300*304*302 

PH2 

1720 

302 

IF( AMX(K) J  9900,304,8800 

PH2 

1730 

8800 

IF(J— 1)9900,303 *8801 

PH2 

1740 

8801 

KP=K-IMAX 

PH2 

1750 

IF(AMX(KP)1 9900* 8803.303 

PH2 

1760 

A  CHECK  E  TO  INSURE  THAT  THE  BOTTOM  ZONES 

Of  THE  PROJECTILE  EMPTY  (FOR  HYPERVELOC ITY)  UP  UNTIL 

8803 

IF( AIXt  K )-Z4 122) )350 ,303,303 

350 

IF(AMX(C) ) 9900 » 303 *306 

303 

M=K 

PH2 

1780 

JJ*J 

PH2 

1790 

GO  TO  307 

PH2 

1800 

304 

AMPY=0.0 

PH2 

1810 

308 

AMUT*0.0 

PH2 

1820 

AMVT=0.0 

PH2 

1830 

0£LET*0.0 

PH2 

1840 

GO  TO  5C1 

PH2 

1850 

300 

IF ( VEL 19901 ,305,304 

PH2 

1860 

305 

If ( AMXILl 19903# 304* 306 

PH2 

1870. 

306 

N=L 

PH2 

1880 

JJ=J*1 

PH2 

1890 

307 

If  (VEU  6130, 6130,6140 

PH2 

1900 

6130 

WSA=(V(Ki*V(L) 1/2.0 

PH2 

1910 

WSB=1.0*(VU.)-V(K))/(DY( JJ ) *SB0UN01 *DT 

PH2 

1920 

VABOVEs WSA/MSB 

PH2 

1930 

CALCULATE  THE  MASS  FLUX  AT  THE  TOP  OF  CELL  K. 

6140 

AMPY=AMX( Ml  * VABOVE/DYi J J ) *0T 

PH2 

1940 

501 

IF(URR) 500*504*502 

PH2 

1950 

502 

If (AMX(K) 19900, 504,503 

PH2 

1960 

503 

IF ( I-ll 7001 , 7001 , 7002 

7001 

M=K 

N=I 

PH2 

1980 

GO  TO  508 

PH2 

1990 

7002 

KP=K— 1 

IFIAMX(KP) 19900, 7003, 7001 

7003 

1F(AIX1KJ-Z( 1221 1357*7001, 7001 

357 

If (I-IMAXJ 351, 7001, 7001 

351 

If (AMX(K+11 19900,7001,507 

504 

AMMP-0.0 

PH2 

2000 

AMUR-0.0 

PH2 

2010 

AMVR=0.0 

PH2 

2020 

0ELER=0.0 

PH2 

2030. 

GO  TO  1 

PH2 

2040 

uuo 
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,  500  IF(FSi9905, 506,504 

506  IF(AMX<K+i) 19904,504,507 
5C7  M=K+1 
N=l  +  1 

5C6  IF £ FS 16100, 6100, 61 10 
6100  WSA=(UIKH-U(K  +  U  1/2.0 

WSB=1.0+I U(K*11-U( Ki l/( DX£N1*SB0UND1*DT 
URR=hSA/WSB 

C  CALCULATE  THE  MASS  FLUX  AT  THE  RIGHT  OF  CELL  K. 

6110  DEN*AMX (Ml/TAUI N) 

IF (GAM1 9989, 9990, 9989 

9989  CART* 1. 

GO  TO  9991 

9990  CART=X( I ) / • 5  / 

9991  AMMP=DEN/PirTS*CART*URR 

1  I F( AMMP1 16 , /  ,8820 

8820  IFIGAMCIJ1 174,74,8821 

8821  1F(FS16120,6120,74 

6120  IFIAMXIK+11 19903,8822,74 

8822  IF(AMX(K1/{TAUU)*DYIJ11-ZI1111  18823,74,74 

8823  1FIAIX<K1-Z( 110118824,74,74 

8824  W S=GAMC 4  J I ♦AMXl K 1-TAUl 1 1 *DY ( J 1 *Z( 111 1 

4  IF ( WS  18  826, 8826, 882  5 

8825  AMMP*HS 
GO  TO  74 

.  8826  AMMP*0. 0 
74  JTAG*0 

BEGIN  CHECKING  TO  SEE  IF  THEIR  IS  ANY 
PREFERENTIAL  MASS  FLUX  BECAUSE  OF  CHOICE  OF 
INDEXING  DIRECTION. 

2  IF( AMPY 13,4,4 

C  TOP  FLUX  IS  INTO  CELL  K. 

3  I  TAG=1 
WSB=AMPY 
AMPY=Q.O 
GO  TO  64 

4  I  TAG=0 

64  IF ( AMMY 19,5,5 

C  BOTTOM  FLUX  IS  INTO  CELL  K. 

5  IF(GAMC(Ji 17,6,6 

C  LEFT  FLUX  IS  INTO  CELLK. 

6  WS=AMXI K 1 
GO  TO  11 

C  LEFT  FLUX  IS  OUT. 

7  WS*AMX( Kl ♦GAMC ( J 1 
GO  TO  11 

C  BOTTOM  FLUX  IS  OUT  OF  CELL  K. 

9  IFTGAMC(J) 110,8,8 
-C  LEFT  FLUX  IS  INTO  CELL  K. 

8  WS=AMX(K)+AMMY 


PH2  205i 
PH2  206( 
PH2  207( 
PH2  208( 
PH2  209( 
PH2  210 
PH2  211 
PH2  2 12c 

PH2  213C 


PH2  2151 
PH2  2 16C 
PH2  217L 
PH2  2181 
PH2  2191 
PH2  220  C 
PH2  22 1C 
PH2  222 C 
PH2  223C 
PH2  224C 
PH2  225C 
PH2  226c 


PH2  22 7C 

PH2  228c 
PH2  229c 
PH2  230 
PH2  231c 
PH2  232 c 
PH2  233 « 

PH2  234 l 

PH2  235' 
PH2  236l 

PH2  237C 
PH2  238C 

PH2  239C 

PH2  240^ 
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GO  TO  11 

C  LEFT  FLUX  IS  OUT  CP  CELL  K. 

10  WS=AMX(K)*GAMC(JJ+AMMY 

11  WSA-AMPY+AMMP 

12  IF4iiSA-WS175,75,13 
CHANGE  TOP  AND  RIGHT  FLUX  TO  BE  THE 
OLD  FLUX  TINES  THE  MASS  OF  THE  CELL/THE  SUM 
OF  THE  OLD  FLUXES. 

13  AMPY»AMPY*WS/bSA 
AMMP»AMMP*WS/WSA 

75  IF( JTAG114, 73,14 
73  WSCsAMMP 

14  IFI 1TAG115« 7000, 15 

15  AMPY*WSB 
ITAG=0 

C  GO  CHECK  CELL  ABOVE. 

GO  TO  40 

C  RIGHT  FLUX  IS  INTO  CELL  K. 

16  IFI FS ) 76, 17*  76 

76  NSC*AMMP 

C  I=IMAX,  SO  CHECK  CELL  ABOVE  K. 

GO  TO  40 

17  I F  C I *l—IMAX 1 19, 18# 9908 

18  URRRSU I K+ 11/2.0 
GO  TO  20 

19  URRR31UIK+1J+UIK+21 1/2.0 

20  IFIURRRI39, 39,21 

C  FLUX  IS  OUT  OF  THE  RIGHT  OF  CELLIK+D . 

21  IFIGAM19992, 9993,9992 

9993  CARTSX( I*ll*2. 

GO  TO  9994 

9992  CART31. 

9994  URRR3URRR/TAUU*114AMXIK+ll/PIDTS*CART 

22  IFI J— 119909,22.24 

23  VBLQ3Vl K4-H /2.0 
GO  TO  26 

24  KP-K+l-IMAK 
VBL03IVlK*U*VlKPll/2.0 

26  IFI VBL0125, 38,38 

C  FLUX  IS  OUT  OF  THE  BOTTOM  OF  CELLIK*11. 

25  VBLO-AMXIK+ll/DYI J1 ♦VBLO*DT 

27  IFI VEL 128, 29, 28 

28  VAB3VIK+11 /2.0 
GG  TO  31 

29  KP3K*IMAX+1 

VAB=t VI K+11+VIKP1 1/2.0 

31  IFI VAB 136, 36,30 

C  FLUX  IS  OUT  OF  TOP. 

30  VAB-AMXIK+ll/DYI J1*VAB*DT 

32  WS-AMXI K4-1 1 


PH2  2410 

» 

PH2  2420 
PH2  2430 
PH2  2440 


PH2  2450 
PH2  2460 
PH2  2470 
PH2  2480  C 
PH2  2490 
PH2  2500 
PH2  2510  c 

PH2  2520 

C 

PH2  2530  c 
PH2  2540 

f H2  2550  c 

PH2  2560  * 

PH2  2570 
PH2  2580 
PH2  259p 
PH2  2600 


PH2  2620 
PH2  2630 
PH2  2640 
PH2  2650 
PH2  2660 
PH2  2670 

PH2  2680 
PH2  2690 
PH2  2700 
PH2  2710 
PH2  2720 
PH2  2730 
PH2  2740 

PH2  2750 
PH2  2760 
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33  WSA=URRR-AMMP— VBLO+VAB 

34  IF(WSA-WS) 77,77,35 

35  AMMP= AMMP*WS/ WS A 
77  JTAG= 1 

KSC=AMMP 
AMMP=Oi 0 
GO  TO  2 

C  FLUX  AT  TOP  IS  INTO  CELL  (K+li. 

36  W  S=AMX ( K+l  I 

3 i  WSA=URRR-AMMP-V8LG 
GO  TO  34 

C  FLUX  IS  IN  FROM  BOTTOM  INTO  CELL  (K*l). 

38  V8LU-=0.0 
GO  TO  27 

C  FLUX  IS  INTO  CELL  IK*1)  FROM  RIGHT. 

39  URRR=C.O 
GD  TO  22 

RIGHT  FLUX  OUT  OF  CELL  (KJ  IS  POSITIVE  AND  TOP 
FLUX  IS  COMING  INTO  CELL  (Ki  FROM  (K+IMAX). 

40  IFIVELI700 0,41, 7000 

41  IF ( FS 142,43 ,42 

C  WE  ARE  AT  THE  RIGHT  BOUNDARY  OF  THE  GRID. 

4  42  KP=K+IMAX 

URT=U( KPI/2.0 
GO  TO  45 

*  43  KP=K+IMAX 

URT*(UIKPJ+U(KF«-1)  J/2.0 

45  I  FI  URT) 46, 46, 70 

C  FLUX  AT  RIGHT  ICELL  Mi  IS  NEGATIVE. 

46  URT =0.0 
GO  TO  47 

70  KP=K*IM4X 

IF (GAMJ 9996, 9997, 9996 

9997  CART=X(I)*2. 

GO  TO  9998 

9996  CART=1. 

9998  URT=URT/TAU( I)*AMX1KPI/PIDTS*CART 

C  FLUX  AT  RIGHT  I CELL  Mj  IS  POSITIVE. 

47  IF( J+1“JMAX)48,49,9910 
4d  KP=K+ IMAX 

KLL=KP+IMAX 

VABT=I V ( KP ) +V (KLL) 1/2. 

GO  TO  51 
49  KP=KHMAX 

KLL=KP+IMAX  / 

VABT=V(KPJ/2.0 
51  IF ( VABT 18810,71,72 

C  FLUX  IS  IN  FROM  TOP  OF  CELL  M.  / 

•  8810  I FI AMX I K) 19903,8811,71 

C  CHECK  FOR  SOLID  OR  VAPOR.  / 


PH2  2770 
PH2  278C 
PH2  279c 
PH2  280i 
PH2  2810 
PH2  2820 
PH2  2830 

PH2  2840 
PH2  2850 
PH2  2860 

PH2  2870 
PH2  2880 

PH2  2890 
PH2  2900 


PH2  2910 
PH2  2920 

PH2  2930 
PH2  2940 
PH2  2950 
PH2  2960 
PH2  2970 
PH2  2980 

PH2  2990 
PH2  3000 
PH2  3010 


PH2  3030 
PH2  3040 


PH2  3070 
PH2  3080 

PH2  3100 
PH2  3110 

PH2  3120 
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8811 

8812 

8813 

8814 

8815 

8816 


8817 

71 

72 

52 

53 


54 

55 


60 

61 

59 

58 

57 

56 

7000 

309 

8833 

8835 

8836 


C 

C 

C 

c 


8837 

8838 


8834 

8839 

8840 

8841 


318 

322 

323 


C 

C 


b 
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if fAMXfKPi/f TAU111*0Y( J+ll )-Z (111)1 88 12, 71, 71 

IF4AIX(KP1-Z(1101I8813,71,71 

VAST®VA8T*AMX (KLL )/DY( J+2) *0T 

WS=~VABT+AMX4KP)-TAUU1*0Y(J*11*Z(1111 

IF(WS 18817, 8817,8816 

AMPY--WS 

GO  TO  71 

AMPY-0.0 

VABT=0.0 

GO  TO  60 

VAflT*VABT*AMX(KPl/DY(J+l)*OT 
if f GAMC ( J+ 11 154, 53 , 53 
WS-AMX4KP1 
GC  TO  55 

MS=AMX4  KP1  -fOAMC  <  J  +  li 
WSA®VABT— AMPY+URT 
GO  TO  57 

IFfGAMCf J+l) 161,61,59 
WS-AMX4  KPJ  *GAMC<  J  +  l ) 

GO  TO  58 
MSsAMX(KP! 

WSA=-AMPY+U«T 

Iff  WSA— 8S1 7000,7000,56 

AMPYsAMPY*WS/l«SA 

GO  TO  7000 

AMMP~USC 

Iff AMP Y 18834, 8831, 8833 
If ( JMAX-319911,318,8835 
KP=K+IMAX 

If 1AMXI KP1 19900,8837,318 

44#*  NOTE  *************************************************** 
ACROSS  FREE  SURFACE,  HOLD  UP  MASS  FLUX 
UNLESS  THIS  MASS  PRODUCES  A  DENSITY  GREATER  THAN  TOZGNE. 
************************************************************ 
IF4 AMPY/f  TAUf 1 1*DY ( J1 l—TQZ ONE) 8838 ,318, 318 
AMPY=0.0 
GO  TO  8831 

Iff J-ll 9911*325,8839 
IF (AMX(K) 19900, 8840, 325 

If 4— AMPY/4  TAU41 )*DYf  J+ll J—TOZONE 18841, 325, 325 

AMPY»0«0 

GO  TO  8831 

OELM=GAMC ( J 1 +AMMY— AMP Y 
Iff VEL19901,324,323 
WS=U(K1**2«-V(K)**2 

MATERIAL  HAS  LEFT  THE  TOP,  TRIGGER  REZONE 

FLAG,  REMOVE  ITS  ENERGY  FROM  ETH( TOTAL  ENERGY  OF  SYSTF.M1. 

ETH-ETH— AMPY*f AIX ( K) +WS/2.Q1 

If fAMPY/fTAUf I )*DY4J11-VT1 324,324,6900 

REZS 1*0 


PH2  3130* 
PH2  3140 

PH2  3160 
PH2  3170 
PH2  3180 
PH2  3190 
PH2  3200 
PH2  3210 
PH2  3220 
PH2  3230 
PH2  3240 
PH2  3250 
PH2  3260 
PH2  3270 
PH2  3280 
PH2  3290 
PH2  3300 
PH2  3310 
PH2  3320 
PH2  3330 
PH2  3340 
PH2  3350 
PH2  3360 
PH2  3370 
PH2  3380’ 
PH2  3390 
PH2  3400 
PH2  3410 
PH2  3420 


PH2  3430 
PH2  3440 
PH2  3450 
PH2  3460 
PH2  3470 
PH?2  3480 
PH2  3490 
PH2  3500 
PH2  3510 
PH2  3520 
PH2  3530 


PH2  3540 
PH2  3560 


n  n  o  r>  o  r> 
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324  AMUT=AMPY*U(K} 

AMVT=AMPY*V(Kj 
GO  TO  326 

325  CONTINUE 

8831  AMUT=AMPY*U(L) 

AMVT-AMPY*VIL) 

DELM=GAMC ( J J— AMPY+AMMY 

326  1F4AMPY)327,328,328 

32  7  OELET=4iX<Li-MUlLI**2+V(L)**?J/2.0 
GO  TO  333 

328  if (AMMYi329, 330,330 

329  D£LET=DELEB 
GO  TO  333 

330  IF (GAMC(J) 1331,332,332 

331  DELET=SIGC l J  ) 

GO  TO  333 

332  DELETE  IX tKi+iUlK)**2+VIKJ**2)/ 2.0 

SUM  UP  RADIAL  MOMENTA  FOR  ALL  FLUXES  EXCEPT 
THE  RIGHI  AND  STORE  IN  SIGMU. 

333  SIGMU=FLEFTI JJ+AMMU-AMUT 

SUM  UP  AXIAL  MOMENTA  FOR  ALL  FLUXES  EXCEPT  THE 
RIGHT  AND  STORE  IN  SIGMV. 

SIGMV=YAMCi JJ+AMMV-AMVT 

SUM  UP  TOTAL  ENERGY  CARRIED  BY  THESE  FLUXES 
EXCEPT  THE  RIGHT  FLUX  AND  STORE  IN  DELEK. 
DELEK^GAMC ( J J  *S IGC ( J1+AMMY*0ELEB— AMPY*DELET 
509  IFUMMPJ8843, 518, 8844 

8844  IF ( I MAX— I 19911,518,8 345 

8845  IF4AMX(K+1T )9900, 8846, 518 

8846  If (AMMP/ITAUII)*DY(JJ J-TGZ0NEI8847,518,518 

8847  ..^MP=0.0 
GO  70  518 

8843  IF( I“1J  991 1,512,8848 

8848  IF<AMX<1<J  >9900,8849,512 

8849  IF (-AMMP/I TAU(I*1J*DY( J) )-T OZONE  1 8 850, 5 12, 512 

8850  AMMP=0.0 
GO  TO  518 

512  OELM=DELM~AMMP+AMX(K» 

513  CONTINUE 

514  CONTINUE 

8828  AMUR=AMMP*U4K*1J 
AMVR=AMMP*V(K+1) 

GO  TQ  525 

518  DELM=DE LM-AMMP+AMX I K ) 

521  CONTINUE 

522  IFIFSJ9905, 524,523 

523  WS=U{K)**2*V(K)**2 
ETH-ETH-AMMP*! AIX( K  J+WS/2-0) 
IF(AMMP/lTAUU)*DYt  JJ  )-VTl 524, 524, 6901 

6901  REZ=1.0 


PH2  357C 
PH2  35 dC 
PH2  3590 
PH2  360o 
PH2  3610 
PH2  3620 
PH2  3630 
PH2  3640 
PH2  365C 
PH2  366C 
PH2  367C 
PH2  3680 
PH2  369C 
PH2  37CG 
PH2  3710 
PH2  3720 
PH2  3730 


PH2  3740 


PH2  3750 


PH2  3760 
PH2  3770 
PH2  3780 
PH2  3790 
PH2  3800 
PH2  3810 
PH 2  3820 
PH2  3830 
PH2  3840 
PH2  3850 
PH2  3860 
PH2  3870 
PH2  3880 
PH2  3890 
PH2  3900 
PH2  39 It 
PH2  3920 
PH2  3930 
PH2  3940 
PH2  3950 
PH2  3960 
PH2  3970 
PH2  3980 

PH2  4000 
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524  AMUR=AMMP*U ( K ) 

AMVR»AMMP*V(K) 

525  SIGMU=S IGMU— AMUR 
SIGMV-SIGMV-AMVR 

526  TIC=0«0 

527  IFtAKMPi 528,529,529 

528  DELER*AiXIK*iH-lUIK*il**2*V(K+l)**21/2.0 
GO  TO  537 

529  IF( AMMY1530, 531, 531 

530  DELER=OELEB 
GO  TO  536 

531  IFIGAMC(J) J 532,533,533 

532  0£LER=SIGC( J1 
GO  TO  536 

533  IF( AMPY J  535 , 535,534 

534  0£LER=DELET 
GO  TO  536 

535  DEL£R=AiX(Kl+4U(Kl**2+V(Kl**21/2.0 

536  TIC«1.0 

537  OELEK=OELEK-AMMP*DELER 

538  IF (TIC 19907 ,539,550 
550  WS=OEL£R 

GO  TO  999 

539  WS-A1XIK1+ (UiKJ**2+V (K1 **21 /2«0 
999  IF(D£1.M1998, 543,540 

998  IF(AMX(KJ*1.E-6*0ELMJ9906,997,997 
997  DELM=0.0 
GO  TO  543 

ENK=T0TAL  ENERGY  OF  CELL  (K)  +  ENERGY  THAT 
HAS  BEEN  AOOEO  AND  LOST. 

540  ENK^AMX (K1 *WS+DEL£K 

BY  CONSERVING  AXIAL  MOMENTA,  CALCULATE  THE  NEW 
AXIAL  VELOCITY  COMPONENT  FOR  CELL  K. 

541  U(K1=4SIGMU+AMX(K1*U(K1 1/OELM 

BY  CONSERVING  RADIAL  MOMENTA,  CALCULATE  THE  NEW 
RADIAL  VELOCITY  COMPONENT  FOR  CELL  K. 

601  VIK1-4  SIGMV+AMX(K)*V(K) 1/DELM 
IF( I— I 1 1603, 6604,6604 

6604  IF1 U(K« 16605,6606,6605 

6605  NRC”i 

6606  1F( VI KJ 16607,6608,6607 

6607  NRC= 1 

6608  IF (AIX1K)1 6609,6610,6609 

6609  NRC= 1 

6610  CONTINUE 

603  WS-U(K1**2*V(K)**2 

BY  CONSERVING  BOTH  TOTAL  ENERGY  AND 
MOMENV CALCULATE  THE  NEW  SPECIFIC 
INTERNAL  ENERGY  FOR  C^LL  K. 

542  AIXIKl=£NK/DELM-WS/2.0 


PH2  4010. 
PH2  4020 
PH2  4030 
PH2  4040 
PH2  4050 
PH2  4060 
PH2  4070 
PH2  4080 
PH2  4090 
PH2  4100 
PH2  4110 
PH2  4120 
PH2  4130 
PH2  4140 
PH2  4150 
PH2  4160 
PH2  4170 
PH2  4180 
PH2  4190 
PH2  4200 
PH2  4210 
PH2  4220 
PH2  4230' 
PH2  4240 
PH2  4250 
PH2  4260- 
PH2  4270 
PH2  4280 


PH2  4290 


PH2  4300 


PH2  4310 
PH2  432 C 
PH2  4330 
PH2  4340 
PH2  4350 
PH2  4360 
PH2  4370 
PH2  4380 
PH2  4390 
PH2  4400 


PH2  4410 


f 
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543 

AMX( KJ=DELM 

PH2 

44 

IF ( AMXI K) ) 9900, 2007, 544 

PH2 

4*-: 

2007 

AIX{ K )=0.0 

PH2 

44 

UI K) =0. 0 

PH2 

44 

V  IK) =0«  0 

PH2 

4^ 

PIKJ=0.0 

THE  RIGHT  VALUES  OF  CELL  < K 1  BECOME  THE  LEFT 

VALUES  OF  CELL  (K+l). 

PH2 

44 

544 

GAMC  ( J) --AMMP 

PH2 

44 

FLEFTC JJ=AMUR 

PH2 

44 

YAMC ( J)=AMVR 

PH2 

4  5 

SIGCI J1=DELER 

THE  TOP  VALUES  OF  CELL! KJ  BECOME  THE 

BOTTOM  VALUES  FOR  CELL  (K+IMAX1. 

PH2 

45 

545 

AMMY=AMPY 

PH2 

41 

AMMU=AMUT 

PH2 

4!: 

AMMV=AMVT 

PH2 

41 

0EL£B=0£LET 

PH2 

41 

546 

K=K+IMAX 

PH2 

41 

LL=K~*  IMAX 

PH2 

41 

IFI UILLJ 16500,6600,6500 

PH2 

41 

6500 

NRI=1 

PH2 

41 

6600 

IF tVILL 116601,6602, 6601 

PH2 

46 

6601 

NRT=  1 

PH2 

46 

6602 

IF (AIX(C)J6611, 54 7,6611 

PH2 

4t 

662  1 

NRT=1 

PH2 

46 

547 

CONTINUE 

PH2 

46 

11=1 1*NRC 

PH2 

4 6 

12= I2+NRT 

PH2 

46 

IFUMAX-1 1)6700,6701,6702 

PH2 

46 

6700 

1 1= imax 

PH2 

46 

6701 

CONTINUE 

PH2 

46 

6702 

IF! J MAX-12) 6800, 6801 ,6802 

PH2 

47 

6800 

I2=JMAX* 

PH2 

47 

6801 

CONTINUE 

PH2 

47 

6802 

GO  TO  548 

PH2 

47 

9901 

NK=300 

PH2 

47 

GO  TO  9999 

PH2 

47 

9900 

NK=302 

PH2 

47 

GO  TO  9999 

PH2 

47 

9903 

NK=305 

PH2 

47 

GO  TO  9999 

PH2 

47 

9904 

NK=506 

PH2 

46 

GO  TO  9999 

PH2 

46 

9905 

NK=500 

PH?. 

46 

GO  TO  9999 

PH2 

46 

9906 

NK=513 

PH2 

46 

GO  TO  9999 

Prl2 

46 

9911 

NK=8833 

PH2 

46 

GO  TO  9999 

PH2 

46 
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9908  NK=  17 

GO  TO  9999 

9909  NK=  22 

GO  TO  9999 

9910  NK=  47 

GO  TO  9999 
9907  NK=538 
9999  NR-4 

CALL  DUMP 
548  SUM=0„0 
2005  DO  2001  I«i,Jl 
K=  I  + 1 

DO  2000  J= 1 , 12 
IFIAMX(K) 12000,2000,2009 

C  If  ANY  RHO  (CELL  DENSITY]  IS  LESS  THAN  TOZONE, 

C  SET  THE  MASS  TO  ZERO,  AND  TALLY  THE 

C  MOMENTAS  AND  ENERGIES  IN  THE  Z  STORAGE,  ALSO 

w  CHECK  FOR  NEGATIVE  INTERNAL  ENERGIES,  IF 

C  WE  FIND  SOME,  SET  THEM  TQ  ZERO  AFTER 

C  SUBTRACTING  THEM  FROM  ETH.. 

2009  IF  I AMX (KJ/ITAUlij *DY ( J ) )-TOZONE 120 10, 2008, 2008 

2010  WS=lU(K]**2+VlK)**2J/2.0 
ZI100)=Z< iOOi+AMXlK) 

WS=AMX(K)*1AIX(K]+WS) 

Z(1Q11=Z(101J+WS 

ETH=ETH-WS 

Z<102)^ZU02)+AMXIKJ*um 

Z I 103  )  =  ZI103) +AMX1K) *V  I K  J 

AMX(K)=0.0 

A I XI K  J  =  0.0 

PIK)=0.0 

UIK)=0.0 

V(K)=0.0 

GO  TQ  2000 

2008  IFIAIXIK) )2004, 2000, 2000 
2004  SUM=SUM+AI X ( K )*AMX( K1 
AIXI K)=0«0 

2000  K=K+IMAX 

2001  CONTINUE 
2003  ETH=E  TH-SU.M 

ZI 104i“ZI 104I+SUM 
800Q  IFIREZJ8001, 8001, 8002 

8002  IFIREZFCT ) 300A, 8004,8003 
8004  REZ=0- 

GO  TQ  800' 

8003  CALL  REZONE 
8001  RETURN 

END 


PH2  4880, 
PH2  4890 
PH2  4900 
PH2  4910 
PH2  4920 
PH2  4930 
PH2  4940 
PH2  4950 
PH2  4960 
PH2  4970 
PH2  4980 
PH2  A 990 
PH2  5000 
PH2  5010 


PH2  5020 
PH2  5030 
PH2  5040' 
PH2  5050 
PH2  5060 
PH2  5070* 
PH2  5080 
PH2  5090 
PH2  5100 
PH2  5110 
PH2  5120 
PH2  5130 
PH2  5140 
PH2  5150 
PH2  5160 
PH2  5170 
PH2  5180 
PH2  5190 
PH2  5200 
PH2  5210 
PH2  5220 


PH2  5260 
PH2  5270 
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* 


1BFTC  £S  LIST ,DcCK»  REF 

SUBROUTINE  ES 


ES 

0900 

METALLIC  EQUATION  OF  STATE,  SEE 

GA-3216  REPORT. 

ES 

0980 

10 

RHOW= AMX(K ) /(TAU(I)*DY(Ji) 

ES 

0990 

LTA=RHOW/Z( 115) 

ES 

1000 

V0W=1. O/ETA 

ES 

1010 

n 

P1=A IX( K) *RHGW*Z1 1  lt> ) 

ES 

1020 

;2 

P2=IZ(115)*IAUI I )*0Y(J)  )**2*A1X(K) 

ES 

1030 

13 

P3=AMXI Ki**2*ZI 117) 

ES 

1040 

14 

P4=Z ( 1 1 8) / ( P2/P3+1 .0 ) *AIX ( K ) *RHOW 

ES 

1050 

15 

P5=Z( 119)*( ETA— 1.0) 

ES 

1060 

16 

IFIET A— 1.0)50,100, 100 

ES 

1070 

50 

IF! VOW— Z I 120) ) 55, 55 ,  .5 

ES 

1080 

55 

I  FI AIX { Ki— Z (122 1)100,100,75 

ES 

1090 

75 

P7=Z(123)<‘(  VOW-l.O) 

ES 

1100 

IF (P7-88. 0)4002, 4002, 4003 

ES 

1110 

4003 

P7=88.0 

ES 

1120 

4002 

CONTINUE 

ES 

.30 

P8=EXP1P7) 

ES 

1140 

P9=l®  0/P8 

ES 

1150 

P 10=Z ( 124) *1 (VOW-1. 0)**2) 

ES 

1160 

IF (P10-88. 0)4000, 4000, 4001 

ES 

1170 

4001 

P 10=88® 0 

ES 

1180 

4000 

CONTINUE 

ES 

1190 

P11=EXPI P10) 

ES 

1200 

P12= 1.0/P1 1 

ES 

1210 

PIK)=P1+(P4+P5*P9)*P12 

GO  TO  119 

ES 

1220 

100 

PS=Z( 126)*I I  ETA-1. 0)**2) 

ES 

1230 

P ( K )=P1 +P4*P5+P6 

ES 

1240 

119 

IF ( PIK) ) 999, 999 5 200 

ES 

1250 

200 

WSGX=. 5 

ES 

1260 

GO  TO  500 

ES 

1270 

999 

P(K)=0.0 

ES 

1280 

WSGX=.5+Z( 125) 

ES 

1290 

*0  TO  500 

ES 

1300 

500 

RETURN 

ES 

1310 

ENO 

ES 

1320 

u  o  o  u  o 
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SIBFTC  REZONE  L I  ST , DECK , REF 
SUBROUTINE  REZONE 

CONSERVE  MOMENTUM  ANO  TOTAL  ENERGY, 
ALL  LINEAR  DIMENSIONS  BY  2.  1THUS  4 
IN  THE  OLD  GRID  ARE  COMBINED  INTO  1 
THE  NEW  GRID.! 

NIMAX=I MAX/2 
NJMAX= JMAX/2 
DO  10  J=l, NJMAX 
K=( J— 1 ) *NIMAX+2 
L=(  J— 1 ) *2* I MAX+2 
DO  11  1  =  1 , NiMAX 
M=L+I MAX 

12  WSA=AMXlL)+AMXIMJ+AMX(L+l)+AMX( M+1I 


AIXI M )=0. 0 

P ( M ) =0  o  0 

AMXIL+i J=0.0 
U ( L+ 1 ) =0*0 
VIL+1J=0„0 
AIX«L+l)=0.0 
P (L+l J=0.0 
AMX(M+1 1=0*0 
U  ( M+l )  =  0  *0 
V(  M+l )=0.0 
AIXIM+11=0.0 
P { M+l J=0«0 


REZOOO 10 
REZ00980 

INCREASE 

CELLS 

FOR 

REZ00990 
REZ01000 
REZ01010 
REZ01020 
REZ01030 
REZ01040 
REZ01050 
REZ01060 
REZ01070 
REZO  1080 
REZ01090 
REZ01100 
REZQ1110 
REZQ1120 
REZ01130 
REZQ1140 
REZO  11. '30 
REZ01160 
REZ01170 
REZ01180 
REZ01190 
REZ01200 

REZ01210 
REZ01220 
REZ01230 
REZ01240 
REZ01250 
REZ01260 
REZO 1270 
REZ01280 
REZO 1290 
REZ01300 
REZ01310 
REZ01320 
REZ01330 
REZ01340 
REZ01350 
P.EZ01360 
REZ01370 
REZ01380 
REZ0139-0 
REZ01400 


WSB=AMX1L)*IUILJ**2  +  VIU**2)+AMX(M)*IU(M) 

1**2+ V(  M 1**2 ) +AMX IL+ 1 1 *( U( L  + 1 i **2+V (L  +  l i **2i 
2+AMX (M+1)*IUIM+1)**2+VIM+1)**2) 

UlK)  =  <UlLJ*AMXlL)+UlMJ*AMXlM)+UlL  +  l)*AMXIL  +  n* 
1U(M+1J*AMX(M+1))/WSA 

V(K)=( VlLl*AMX(Li+V(MJ*AMX(M)+V(L+ll*AMX(L+i)+ 
1V(M+1)*AMXIM+U )/WSA 

AIX ( K J  =  A I/I L J  *AMX<  L  J  +  A1X l M) *AMX  (M)+AIX(L+1)* 
lAMXlL+l)+AMX(M+l)*AIXiM+l) 

AMXlK)  =  i.SA 
WS=U(KJ**2+V(KJ**2 
E=A1X(K L+WSB/2.0 
AIXIKI=E/AMX(KI— »5*WS 
IFIK-2J 14,14,13 

C  SET  CELL  QUANTITIES  OF  OLD  GRID  TO  ZERO. 

13  AMXIL)=0.0 
U(L)=0*0 
V (Li =0. 0 
AIX(L )=0.0 
P ( L I =0.0 
AMX ( M )=0.0 
U ( M  )  =  Q*  0 
V(M)=0«Q 


14  K=  K  + 1 
L=L+2 

11  CONTINUE 
10  CONTINUE 

C  CALCULATE  NEK  DY  ANU  Y  I  UMAX  OF  THEM]. 

18  DO  999  J= 1 ,  JMAX 
DY( J)=OY( J)*2.0 
999  CONTINUE 

DO  99  J  =  1 »  JMAX 
YiJ)=Y( J-1)+DY( J) 

99  CONTINUE 

16  DX(1)=2.0*DX11) 

XI  1)=DX(1) 

KS=X( l)**c 

IF (GAM) 3001,3000, 3001 

3001  KS=DX(1) 

3000  TAU( 1)=PIDY*KS 

C  CALCULATE  NEK  DX  AND  X,  AND  TAUIIMAX  OF  THEM! 

17  DO  98  1=2, IMAX 
XII)=XI I— 1 )  *DX( 1 ) 

Dxm=oxm 

KSA=X1I )**2 

IF (GAM) 3002, 3003, 3002 

3002  TAU(I)=DX(I) 

GO  TO  98 

3003  TAU< IJ=PIDY*(KSA-WS) 

WS=WSA 

98  CONTINUE 
IMAX=NI MAX 
JMAX=NJMAX 

C  PREPARE  NOW  TO  SHUFFLE  THE  K  ARRAYS  SUCH 

C  AS  TO  PRESERVE  K=<  J- U*  I  MAX*  I  ♦  1 . 

DO  20  N=1  * JMAX 

J= JMAX+1— N 

K=I  J“1J*IMAX+1-HMAX 

L  =  ( J— 1 ) * ( IMAX  +  IMAXUl  +  IMAX 

DO  21  I  =  i , I  MAX 

1000  AMXI L ) =AMX( K ) 

A 1X1  L  )  =  Al  X  (  K  ) 

U(L)=U(K) 

V ( L)=V( K) 

P(L)=P(K) 

IF  I J— 13 1002,1002,1001 

1001  AMX ( K ) =0.0 
A  IX ( K )  =  0.0 
V(K)=C.C 

U ( K ) =0. 0 
P I K)  =  0«  0 

1002  K=K-1 
L  =  L~ 1 


REZO 14 
REZ014 
REZQ14 
REZ014 

REZ014 

REZ014 

REZ014 

REZ014C 

REZ014' 

REZOi 5v 

REZ015 

REZO 15. 

REZ015 


REZ015 
REZ015 
REZO  15 
REZ015 


REZ0l6t 

REZ016 

REZ016. 

REZOlfc 


REZ016 

REZ016 

REZOlo 

REZ016 

REZ016 

REZOlfa 

REZ017 

REZ017 

REZ017 

REZ017 

REZ017 

REZ017 

REZC17< 

REZ017 

REZ017 

REZ017 

REZ018 

REZ018 
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21  CONTINUE 
20  CONTINUE 

IMAX=NIMAX*2 
JMAX=NJMAX*2 
I I=NI MAX+1 


RE201820. 
REZ01830 
REZO  1840 
REZ01850 


JJ=NJMAX*1 

C  A00  ON  NEW  NASS  WITH  OENSI TY=Z ( 111 )  IN  TAR6ET 


00  50  i=l«NIMAX 
K=( JJ— 1 )*iMAX+I+i 
00  60  J=JJ  « JMAX 
AMX<K)=Z< 111)*TAU<I)*0Y(  Ji 

60  K=K+ IMAX 
50  CONTINUE 

JJ=(Z(147)/2*+*2) 

JJ=JJ*1 

00  61  I  —  1 1 • IMAX 
X*I+l+( JJ-1I*IMAX 
00  62  J= JJ •  JMAX 
AMX<KJ=Z<111)*TAUIIJ*0Y< Jl 
62  K=K*IMAX 

61  CONTINUE 

RESET  ACTIVE  GRIO  MARKERS* 

ASSUMPTION  THAT  ALL  OX  AND  OY  ARE  = 
NOTE*  DVK=K(0)/{RH0I0)*DX  SQ.  i 

ws=oxm*oxm 

OVK=DOXN/(Z( llli*WS) 

C 

JJ*JJ-1 
Z ( 147J=JJ 
Il=NIMAX*2 
I2=N JMAX+2 
WS=T*OTNA 
NK=NC+1 

C  EOIT  THE  NEW  GRID* 

WRITE  < 6, 8004 ) WS »NK, DX( 11 

WRITE  1 6* 8007 UMAX* (X(1),I=0, IMAX) 

WRITE  16*8008) JMAX* 1 Y  ( J  i ,J=0,JMAX) 

WRITE  ( 6,80091 IMAX, I OX<I)» 1=1, IMAX) 

WRITE  (6, 80 101 JMAX, ( DY( J ) , J=1 , JMAX ) 

WRITE  <6,801 l) IMAX, I TAUII), 1=1, IMAX) 

KMAX- I MAX* JMAX + 1 

IMAXA= I MAX+1 

JMAXA= JMAX+ 1 

KMAXA=KMAX+1 

RETURN 

80040FORMATUH  ////22H  PROBLEM  REZONED  AT 
1X-E12.6////) 

8007  FORMAT < 1H  /10H  XII)  1=0 • 12/ < 5F16.6 ) i 

8008  FORMAT ( 1H  / 10H  Y<J)  J=0 , 12/ < 5F16.6 ) I 

8009  FORMAT! 1H  /11H  DX<li  1=1 , 12/ < 5F16. 6) ) 


REZ02040 

REZQ2050 

REZ02060 

REZ02070 

REZ02080 

RE202090 

REZ02100 

REZ02110 

REZ02120 

REZ02130 

REZQ2140 

REZ02150 

REZ02160 

REZ02170 

REZ02180 

= 1PE12*6, 6X ,5HCYCLEI4,6X, 3H0REZ021 90 

REZ02200 

REZ02210 

REZ02220- 

REZ02230 
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8010  FORMAT  t  1H  /11H  OYiJ)  J= 1 • 1 2/ ( 5F 16.6) ) 

8011  FORMAT (  1H  /13H  AREA ( 1 )  1  =  1 1 12/ { Fl6.6t4F  18.6)  )# 
ENO 


REZ0224 

REZ0226 


o  o  o  nr nnn  n  o  n  o  o  o  o  o 
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$ IBFTC  EDIT  LIST, DECK, REF 

SUBROUTINE  EDIT 


SENSE  LITE  (1)  INDICATES  LAST  CYCLE  OF  THIS 
RUN* 

SENSE  LITE  (3)  INDICATES  FIRST  CYCLE  OF  THIS 
RUN. 

104  CALL  SL ITET 1 3»K000FX) 

GO  TOI 106,106) ,K000FX 
106  CALL  SLITE  (31 
GO  TO  126 

108  IF (CYCLE-CST0P1 110,122,122 
110  IF(REZ)9901, 112,124 
112  IFIAHOD (CYCLE, DUMPT 71 1114,124, 114 
114  IF( AHOD (CYCLE, PRINTL 11120,126,120 
120  IF( AH0D(CYCLE«PRINTS1 1 140,128,140 
NORMAL  STOP  UN  THIS  CYCLE. 

122  CALL  SLITE  (1) 

DUMP  ON  TAPE  7. 

124  GO  TO  1 

126  CALL  SLITE  (41 

128  GO  TO  6000 

130  GO  TO  1000 

132  CALL  SLITET (4«K000FX! 

GO  T0( 134,136) ,KOOOFX 
134  GO  TO  5000 

CHECK  FOR  ENERGY  CHECK  ERROR.  WHERE 
ECK=  PERCENT  ERROR/PER  CYCLE. 

ECK*( ETH-EJ/ETH  AT  CYCLE  N-( ETH-El/ETH 
AT  CYCLE  N-NPC  ALL  DIVIDED  BY  NPC,  NOTE 
NPC-  NO.  OF  CYCLES  BETWEEN  ENERGY  CHECK 
136  1F(ABS( ECK1— DMIN1 140,140,9905 
140  CALL  SLITET ( 1»K000FX) 

GO  TOI 142, 144) .KOOOFX 
142  REWIND  7 

CALL  SLITE  (1) 

144  GO  TO  10000 


DUMP  ON  TAPE  7 
1  IF( DUMPT 7 130, 3,3 
3  BACKSPACE  7 
WS-555.0 

WRITE  (  71 WS , CYCLE ,N3 
WRITE  (  7) ( Z ( L ) » L=1 ,MZ) 

6  WRITE  (  7) ( U( K) ,V(K) , AMX( K) ,A1X(K),P(K) , K=i ,KMAXA  ) 

7  WRITE  (  71K(0) , (X1K) ,TAU(K) , K=1,IMAX) 

WRITE  (  7) 1 Y( K) ,K=0, JMAX) 

WS-=666.0 


EDIT0010 
EDIT0990 
EDIT  1000 


EDIT  1040 
EDIT  1050 
EDIT  1060 
EDIT 1070 
EDIT 1080 

EDIT  1 100 

EC  IT 1150 

EDIT  1160 

EDIT1170 
EDIT  1180 
EDIT  1190' 
EDIT 1200 
EDIT1210 
EDIT  1220- 
EDIT  1230 


EDI T 1240 
EDIT  1250 
EDIT  1260 

EDIT  1280 
EDIT 1290  C 

EDI T1300  c 

EDIT1310 
EDIT  1320 
EDIT  1330 

EDIT1360  C 


EDIT  1480  c 


n  o 
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WRITE  I  7) WS, WS, WS 
WRITE  (6»8120)NC 
30  GO  TO  126 


6000  NK=12 

C  TABS  ARE  TANGENT  ALPHAS. 

TAB( 1 )=0.02 
TAB(2)=0.04 
TAB( 31=0.06 
TA8(4)=0.08 
TAB(5)=0.10 
TAB( 6 )=0. 15 
TAB! 7J=0.20 
TAB( 8)=0.25 
TAB(9)=0.30 
TAB( 101=0.4 
TAB( 11)=0.5 
TAB( 12) =1.0 
6010  00  6012  1=1 » 1 8 
6012  PR( I ) =0.0 
NKl=NK+2 
00  6014  1=1, NK1 
AMK( I )=0.0 
PK( I ) =0  .0 
6014  QK(I)=0.0 

00  6028  K=2,KMAX 
6017  PR( 1) =0.0 
PR(2)=0.0 
PR(4)=0. 

C  CALCULATE  KINETIC  ENERGY  IN  CELL  K. 
WSB=(U(K)**2*VIK)**2)*.5 

6019  IF(AMX(K))9917, 6028, 6020 

6020  I=MK1 

IF(V(K) 16026,6026,6022 
6022  WSA=ABS(U(K) )/V(K) 

DO  6024  1=1, NK 

C  SEARCH  FOR  TAN  ANGLE  THAT  VELOCITY  VECTORS 

C  MAKE. 

IFiTABI ii— WS A >6024, 6026, 6026 
6024  CONTINUE 
I=NK+1 

6026  W S=AMX( K) 

C  SUM  UP  MASS  BETWEEN  ANGLES. 

6027  AMK( I )=AMK ( 1 )+AMX(K) 

C  SUM  UP  RAOIAL  MOMENTA  IN  THE  ANGLES. 

PK(I)  =  PK(I)4-UIK)4AMX(K) 

C  SUM  UP  AXIAL  MOMENTA  IN  THE  ANG;.  .  . 

QK( I ) =QK( I ) * V ( K ) *AMX ( Kj 
C  SUM  UP  TOTAL  INTERNA  ENERGY 


EDIT1500 
EDIT151C 
EDIT 152C 
EDIT153C 
EDIT  1540 

EOIT1550 
E0IT1560 
EDIT 157C 
EDIT1580 
EOIT1590 
EO IT  1600 
ED1T1610 
EDIT1620 
EOIT1630 
EDIT  1640 
EDIT1650 
EDIT  1660 
EDIT  1670 
E0IT1680 
EDIT  1690 
EDIT1700 
EDIT1710 
EDIT1720 
EDIT1730 
EDIT1740 
EDIT  1760 
EDIT1770 


EDIT1790 
EDIT  1800 
EDIT  1810 
EDIT1820 
EDIT1830 


ED1T1840 
E0IT1850 
EDIT 1860 
EDIT  1870 

EDIT 1880 

EDIT189C 

EDIT  1900 


o  n 


EDIT  1910 


PR( 5)=PR{ 5) +AIX( K) *AMX( K) 

C  SUM  UP  TGT  A.L  KINETIC  ENERGY 
PR(6)=PR( 6) +WSB*AMX( Ki 
C  SUM  UP  TOTAL  MASS 
PR48)=PR48)+AMX(X) 

6028  CONTINUE 

PRO)  =PR  ( 1 )  *PR  ( 2 ) 

PR47)=PR15)*PR(6) 

XNRG=PR4  7) 

PR(9)=PR(1)*PR(5‘ 

PR I 10)=PR(2)+PR(6) 

PRI11)=PRI3)+PRI7) 

PR4 12)=PR( 4)  +PR ( 8J 
WSA=( ETH-PR4 11) l/ETH 
IF (CYCLE) 9931(9931 ,9932 

9931  NPC=1 

9932  PR 4 18)=(WSA— DNN) /FLOAT (NPC) 

ECK=PR( 18) 

DNN=WSA 

C  RESET  CYCLE  COUNTER  BETWEEN  ENERGY  CHECK. 

NPC=0 
SUM=0.0 
00  800  I=l(13 
SUM=SUM*QK( I ) 

800  CONTINUE 

C  RA0ET=  TOTAL  POSITIVE  AXIAL  MOMENTUM  IN  GRIO 
RADET=SUM 

801  SUM=0.0 

00  810  K=2(KMAX 
IF(AMJUK))  810,810,802 

802  IF(U( K) 1810,810,803 

803  SUM*$UM*AMX(K , ♦UIK) 

810  CONTINUE 

C  RADER3  TOTAL  POSITIVE  RA01AL  MOMENTUM  IN  GRIO. 
RAOER=SUM 
SUM=0.0 
JJ=Z( 147) 

00  811  1=1, I MAX 
K=I*1 

00  813  J=1,JJ 
IF( AMX( K ))813,813,814 
814  IF(U(K) 1813(813(816 
816  SUM=SUM*U4K) *AMX4K) 

813  K=K*IMAX 

811  CONTINUE 

C  RAOEB=  TOTAL  POSITIVE  RADIAL  MOMENTUM  BELOW 
C  INITIAL  TARGET-PROJECT ILE  INTERFACE. 

RADEB=SUM 
PR( 191=0.0 
DO  6029  1=1, NK 


EDIT  1920 

EDIT) 930 
EO IT  1940 
E0IT1950 
EDIT I960 
EO I 71970 
EDIT  1980 
EDIT1990 
EDIT  2000 
EDIT2010 
EDIT2020 


E0IT2040 

EDIT2050 

EDIT 2060 

E0IT2070 

EDIT2080* 

E0IT2090 

EDIT2100 

EDIT2110 

E0IT2120 

EOIT2130 

E0IT2140 

E0IT2150 

E0IT2160 

E0IT2170 

E0IT2160 

E0IT2190 

E0IT2200 

E0IT2210 

E0IT2220 

E01T2230 

EOIT2240 

EDIT2250 

E0IT2260 

E0IT2270 

E0IT2280 


E0IT2290 
E0IT2300- 
EO IT  2310 


r>  o  o 
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6029  PR( I+19)=PR(I+18)+AMK(  I) 

PRINK +20)=0.0 
PRI NK+21)=0.0 

WRfTE  16,8116)PRQB,NC»T,DTNA, TRAD*DTRAD»NR»N1»N2» 
WRITE  I6,8117HPRI i) *  1  =  1*8) 

WRITE  I6,8118HPR(I),I=9,12 

WRITE  (6,8119)RAD£B,RAD£R,RAD£T,UVMAX,ETH,ECK 
WRITE  ( 6, 90401 N 10, Nil, II, 12,13,14 
WRITE  16, 8124) I  I  ,AMKIU,PR(  14-19 )  ,PKi I ) «  OKI  I  ),I  =  1, 
6090  GO  TO  130 

£***«  £MD  Qf  s  P  SUBROUTINE  ************************ 


****  SUBROUTINE  PLOT  ********************************* 


E0IT232 
EDIT233! 
EDIT234 ( 

N3,N4  E0IT2351 

EDIT  236 
EDIT237( 
E0IT2381 
EDIT239i 

NKl )  EDIT240L 

ED1T241G 

*#***#*#*********£DIT242< 

EDIT2431 

E0IT244( 

*♦♦*#♦#******#*♦*£01X2451 


1000 

GO  TO  1030 

EDIT246( 

1030 

WRITE  I 6, 81 16 ) PROB*NC, T , OTNA, TRAO ,0TRAD,NR,N1 ,N2, N3, N4 

EDIT247C 

JMAX= JMAX 

EDIT248t 

WRITE  (6,8307)X1,X2,XMAX,Y1«Y2«Y( JMAX I 

EDIT2491 

M=1 

EDIT250C 

IF! JMAX-52 ) 1034 ,1036,1036 

EDIT251 1 

1034 

M=I ABSI 51— JMAX I /2 

EDIT252( 

1036 

DO  1040  1=1, M 

EDIT253C 

WRITE  16,8308) 

£  D I T 2  54 1 

1040 

CONTINUE 

EDIT2551 

1044 

J=  1 2 

ECIT2561 

1100 

K=( J-11*IMAX+1 

EDIT257t 

1105 

DO  1180  1=1,11 

EDIT2581 

K=K  +  1 

EDIT259. 

C 

REPLACE  600000000000  BY-17179869134 

EDIT260( 

1126 

PR UJ=(-ABSI-17 179869184)) 

EDIT261 ( 

1150 

IF  (AMXIK) 19917, 1166, 1160 

EDIT2621 

C 

X  PARTICLE  ONLY 

EDIT263C 

C 

REPLACE  6700000000  BY  922746880 

EDIT  264 1 

1160 

PRII )=0RIPR4 I ),  ABSI  922746B80)  ) 

ED1T2651 

GO  TO  1180 

EDIT  266 

C 

REPLACE  6000000000  BY  805306368 

EDIT  267  ( 

1166 

PRU)=GR(PRU ),  ABSI  805306368)  ) 

EDIT266 

1180 

CONTINUE 

EDIT  269 

1200 

IF (MODI J,5)) 1210,1204,1210 

EDIT270 

1204 

IF (DYIJ)-OY(J-l) 11206,1208, 1206 

EDIT271 

1206 

WRITE  (6,8211)DY1J) ,J,(PR( I ),1=1,I1) 

E0IT272 

GO  TO  1224 

EDIT273 

1208 

WRITE  I6,8201)J,4PR(I ), 1=1, III 

EDIT274 

GO  TO  1224 

EDIT  275 

1210 

IF (DY(J)-GYI J-l)) 1212, 1214, 1212 

EDIT276 

1212 

WRITE  (6,8222)DYIJ)*(PRII)»I“1»I1) 

EDIT277 

GO  TO  1224 

EDIT278 

1214 

WRITE  (6, 8202J(PR(I) ,1=1,11) 

EDIT279 

1224 

J=J-1 

EDIT280 

1226 

IF!  JH230, 1230,1100 

EDIT281 

88 


1230 


1240 

C**** 

C 

C 

C**«* 

5000 

5004 


REPLACE  604000000000  BY- 17716740096 
PR( 1)=(-ABS1 -17716740096) ) 

WRITE  (6, 8201JJ,(PRU),  1=1,111 
WRITE  16, 8302H 1,1=0, IMAX, 5) 

GO  TO  132 


ED1T2820. 

E01T2830 

E0IT2840 

ED1T2850 

ED1T2860 


END  OF  PLOT  SUBROUTINE  **************♦****♦♦********♦***********£0112870 

EDI T2880 
EDIT2890 

SUBROUTINE  L  P  ***************************************************EDIT2900 


WRITE  16,81 16 )PROB ,NC, T , DTNA,TRAO»  DTRAD,NR, Nl ,N2* N3,N4 
DO  5050  1*1,11 
CALL  SLITE  14) 

J=I2*1 

K=I2*IMAX+1+I 
DO  5046  L=1,I2 
J=J-i 
K=K— I MAX 

5012  IF ( AMX( K) ) 9917,5046,5014 
5014  CALL  SLITET 14, KOQOFX) 

GO  TO! 5016, 5019), KOOOFX 
5016  WRITE  (6, 8135)1, XIII, DX (I) 

C  WS=  DENSITY  OF  CELL! XI  IN  GRAMS/CM.  CUBED. 

5019  WS=AMX(K)/1TAU1I)*DY( J) ) 

C  WSA=  COMPRESSION  *  RHO/RHO  NOT. 

WSA=WS/21 111) 

C  WSC=PRESSURE 

WSC- P ( K J 

C  FIRST  COLUMN*  1J)  THE  ROW  NO. 

C  SECOND  COLUMN  =  RADIAL  VELOCITY  =  U  CM. /SEC 
C  OR  CM./SH. 

C  THIRD  COLUMN  =  AXIAL  VELOCITY  =  V  CM. /SEC 
C  OR  CM./SH. 

C  FOURTH  COLUMN  =  PRESSURE  =  F/A  ERGS/ (CM.  CUBED1 

C  OR  JERKS/ ( CM. CUBED) 

C  FIFTH  COLUMN  =  AMX  =  MASS  IN  GRAMS. 

C  SIXTH  COLUMN  =  RHO  =  OENSITY  IN  GRAMS/CC. 

C  SEVENTH  COLUMN=SPECI FIC  INTERNAL  ENERGY  IN  ERGS/GRAM 

C  OR  JERKS/GRAM 

C  EIGHT  COLUMN  *  COMPRESSION  =  RHO/RHO  NOT 

C  NINTH  COLUMN  =  Z  VALUE  (CM.)  OF  TOP  OF  CELL) 

50180WRI T£  1 6,81081 J ,U(K) ,V(K) ,WSC,AMX(K), 

1WSA, Y1 J) 

5046  CONTINUE 
5050  CONTINUE 
GO  TO  136 


EDIT2910 

EDIT2920 

EDIT2930 

EDIT2940 

EDIT2950 

EDIT2960 

EDIT2970 

EDIT2980, 

EDIT2990 

EDIT3000 

EDIT3010 

EDIT  3020 

EDIT3030” 

EDIT3040 


WS, AIXiK) , EDIT 3060 
EDIT3070 
EDIT3080 
EDIT3C90 
EDIT31O0 


C**** 

C 

c 

c 

9901 


END  OF  L  P  SUBROUTINE  *****************************************£01 T31 10 

EDIT3120 

EDIT3130 

ERROR  EDIT3140- 

NK=1 10  EDIT3150 
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GO  TO  9999 

C  ENERGY  CHECK 

9905  NK= 136 

GO  TO  9999 

C  NEGATIVE  MASS 

9917  NK=6015 

GO  TO  9999 

9920  NK=904 

GO  TO  9999 

9921  NK=912 

GO  TO  9999 

9922  NK=918 

GO  TO  9999 

9923  NK=92 2 

GO  TO  9999 

9924  NK=926 
9999  NR=6 

CALL  DUMP 
10000  RETURN 
C 

C  FORMATS 


ED1T3L6 

E0IT317 

EDIT318 

EDIT319 

EDIT320 

EDIT321 

EDIT322 

ED1T323 

EDIT324 

E0IT325 

EDIT326 

EDIT327 

E0IT328 

EDIT329 

EDIT330 

tDIT33i 

EDIT332 

E0IT333 

ED1T334 

EDIT335 

ED1T336 


8108  F0RMATU3, IX, 1P2E14.6.3E15.6, E14.6.E15. 6,614.6)  E0IT337 

‘  8116 OF OR MAT ( 8H1PRQBLEM6X , 5HCYCLE9X  » 4HT I ME 13X#2HDT 13X, 4HTRA01 IX  »  5H3TR A01EDIT338 
12X,2HNR6X,2HN14X,2HN24Xf2HN34X,2HN4/*F7.i,Ill,2X, 1P4E 16 . 7, 1 10, 2X, 4ED1T339 
21611  E0IT340 

*  81 170F ORMATC 1H0//17X2HA116X, 2HAK 1 4X , 5HA I +AK 1 5X , 2HAM/ 4H  00T3X , 1P4E18. 7/3E0IT341 
1H  X4X,4E18.  7 )  ED1T342 

81180F0RMAT(12X,13H - - - 5X,  13H - 5X,  13H - 5EDIT  343 

IX ,  13H- — - /7H  T0TALS1P4E18.7J  EDIT344 

81190FORMAT ( 2H0  //16X,5HRA0EB13X, 5HRA0ER13X, 5HRA0ET12X ,7HMAX  VEL 13X , 3HTED I T345 
1HE12X, 9HREL  ERR0R/7X  ,  1P6E18., 7////  J  E0IT346 

8120  FORMAT ( 1HQ//2 1H  TAPE  7  DUMP  ON  CYCLEI5////)  EDIT347 

81240FQRMAT13H  K12X.5HAM1 KJllX, 9HSUM  AM(Ki 1 1X,4HP(K) 13X,4HQCK)/( I3,4X, ED1T348 
11P4E18.71  J  EDIT349 

8131  FORMAT* 1H  /// 11H  DY(J)  J=1 , 12// ( 1QF 12. 3 ) )  EDIT350 

8133  FORMAT < 1H  ///11H  Y(J)  J=0, 12/ / ilOF 12. 3 1 1  ED1T351 

81350F0RMAT* 1H  ///4H  1  =I3,6X,6HXU)  =F12. 3 , 6X , 7HDX (I J  =F12.3//3H  J8X,ED1T352 

11HX13X, 1HY13X »3Hf /A12X, 3HAMX12X, 3HRH01 1X*3HA1X12X ,4HC0MP 1 IX, 2H  Y/1EDIT353 

8201  FORMAT ( 1 10 »  2H  I54A2)  EDIT354 

8202  FORMAT' 10X,2H  154A2)  EDIT355 

8211  FCkMATiF7.1,I3,2H  I54A21  EDIT356 

8222  FORMAT* F7ol,3X, 2H  I54A21  EDIT357 

8302  F0RMATU12, 10110)  EDIT358 

83070F0RMAT ( 5H  XI  =  1PE J  2 . 6 »3X , 4HX2  =E12.6,3X,6HXMAX  =E 12 .6 ,6 X , 4HY1  =E12ED1T359 

i.6»3X ,4HY2  =£12.6, 3X.6HYMAX  =E12.6i  EDIT360 

8308  FGRMATUH  /)  ED1T361 

9040  FORMAT* 1H  /  6161  E0IT362 


END 


EDIT363 


rnnnoficnoncrinnpnnoonono^nnnornnonfionnnnonnnnrnnn 
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$ I8FTC  PH3  .,S  DECK , REF 

SUBROUTINE  PH 

NOTE,  THIS  C0^_  IS  CALLED  RPM, 

RIGID  PLASTIC  MATERIAL. 

***  NOTE,  BOUNDARY  CONDITIONS  *** 

AT  FREE  SURFACE  (GRID  BOUNDARY!  SET 
THE  STRESS  AT  THE  RIGHT,  TOP  OR  BOTTOM 
=  TJ  THE  STRESS  AT  THE  LEFT,  BOTTOM  OR 
TOP  RESPECTIVELY,  INSURING  THAT  THE 
ACCELERATION  OF  THE  BORDER  CELLS  BE 
ZERO. 

FOR  AXIS  OF  SYMMETRY,  SET  THE  NORMAL 
AND  SHEAR  STRESSES  TO  0.  SINCE  REGARDLESS 
OF  MANNER  OF  COMPUTING  THEM,  THE  AREA 
OVER  WHICH  THEY  ACT  IS  ZERO. 

FOR  REFLECTIVE  BOUNDARIES  AT  THE  BOTTOM 
OR  IN  THE  GRID,  USING  THE  BOTTOM  AS 
A  EXAMPLE,  SET  U(BOTTOM)  =  UIKi  BUT 
SNT=NORMAL  STRESS  AT  THE  TOP  OF 
A  CELL.  ITHE  ARRAY-ASNI ill. 

ST  .'-SHEAR  STRESS  AT  .HE  TOP  OF 
A  CELL.  ITHE  ARRAY  =  AST I  I ) I • 

STR= SHEAR  STRESS  AT  THE  RIGHT  OF 
A  CEL'..  (THE  ARRAY  =  RST (  1+1 )  )  . 

S NR = NORMAL  STRESS  AT  THE  RIGHT  OF 
A  CELL.  (ThE  ARRAY-RSNI 1+1! ) • 

SNB= NORMAL  STRESS  AT  THE  BOTTOM  OF 
A  CELL.  (THE  ARRAY=ASNBI 1)1. 

STB= SHEAR  STRESS  AT  THE  BOTTOM  OF 
A  CELL.  (THE  ARRAY=ASTB ( I J  3 . 

SNL=NQRMAL  S CRESS  AT  THE  LEFT  OF 
A  CELL.  (THE  ARRAY=RSN( III. 

STL=SHEAR  STRESS  AT  THE  LEFT  OF 
A  CELL.  (THE  ARRAY=P.ST (  1 1 )  • 

ODVK=HOOP  STRESS  OF  A  CELL  ( ARRAY=SIG33 (III. 

X l=UOOT  OF  A  CELL,  1 ARRAY=DUDOT (I  I 

X2=V DOT  OF  A  CELL,  ( ARRAY=DVDOT ( I I 
SIGC(I)=U(K  BELOW! 

GAMC ( 1 ) =V( K  BELOW) 

FEF  IS  A  FLAG  TO  BYPASS  THE  STRENGTH 


PH3  0000 
PH3  0C05 
PH3  0010 
PH3  0015 
PH3  0020 
PH3  0025 
PH3  0030 
PH 3  0035 
PH3  0040 
PH3  0045 
PH3  0050 
PH3  0055 
PH3  0060 
PH3  0065 
PH3  0070 
PH3  0075 
PH3  0085 
PH3  0090 
PH3  0095 
PH3  0100 
PH3  0105 
PH3  0110 
PH3  0115 
PH3  012a 
PH3  0125 
PH3  0130 
PH3  0135 
PH3  0140 
PH3  0145 
PH3  0150 
PH3  0155 
PH3  0160 
PH3  0165 
PH3  0170 
PH3  0175 
PH3  0180 
PH3  0185 
PH3  0190 
PH3  0195 
PH3  0200 
PH3  0205 
PH3  0210 
PH3  0215 
PH3  0220 
PH3  0225 


PH3  0240 


noofio  oono  oonooono  oooo  o  ooo 
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C  CALCULATIONS  PH3  024 

VSAVE=VT 

1FIFEFJ7778,7779,7778 
C  NOTE,  CLEAR  THE  PRESSURE  ARRAY. 

7/79  DO  96G0  K=2,KMAX 
P«K)=0e 
9600  CONTINUE 
WS=  I  3 
DTT=DT 
DT=DT/taS 

C  BIG-DV/OZ  CKIT 

BIG=DVK  *DT*TABLM 

C  PROVISIONS  FOR  SUBCYCLING  THE  STRENGTH 

00  9000  NN= 1 , I 3 

PH3  023 

NOTE,  SET  THE  FLAGS  FOR  INCREASING  THE 
ACTIVE  GRID  COUNTERS  TO  0. 

NRC-=0 
NRT=0 
SUM=0. 


1  DO  3302  J=l, 12 

NOTE,  SET  UP  DO  LOOP  IN  J  DIRECTION  FIRST,  NOTE 
THAT  12  IS  THE  LIMIT,  NOT  JMAX 

2  K=t J-li*IMAX+2 
N=K+ 1  MAX 

NOTE,  K  IS  THE  INDEX  OF  CELL  IN  QUESTION,  N 
IS  THE  INDEX  OF  CELL  ABOVE. 

***,  NOTE,  THE  STRESSES  AT  THE  AXIS  OF  SYMMETRY 
ARE  SET  TO  0.  SINCE  REGARDLESS  OF  MANNER  OF 
CALCULATION,  THE  AREA  IS  0. 

3  SNL=0 «, 

STL=0« 

SET  UP  00  LOOP  IN  THE  I  DIRECTION,  NOTE  THAT 
II  IS  THE  LIMIT,  NUT  IMAX. 

4  DO  3361  1=1,11 


PH3  026 
PH3  026 
PH3  027 
PH3  027 
PH3  028 
PH3  028 
PH3  029 
PH3  029 
PH3  030 
PH3  030 
PH3  031 
PH3  031 
PH3  032 
PH3  032 
PH3  033 
PH3  033 
PH3  034 
PH3  034 
PH3  03b 
PH3  03b 
PH3  036 
PH3  036 


AMDM=  FACTOR  FOR  STRESS  CUTOFF  ON  THE  BASIS  OF 

THE  DENSITY  BEING  LESS  THAN  AMDM  X  THE  INITIAL  DENSITY 


CHECK  FOR  RAREFIED  MATERIAL  IN  CELL  (K) 

40  IF  (  AMX<  K )/(  TAU(  I)*DYI J)  )-AMDM*Z(  115)  )3340,  3340,38 
38  IF! J- 11720, 720, 721 


oor>o  on  nn  o  o  o  o  o  o  o  n  o  n  n 
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721  K BLO=  K-IMAX 

IFlAMX(KBLa)/<TAUU)*DY<  J-U )-AMQM*Z( 11 5{ I  722. 722 ,720 

722  DDVK=0. 

60  TO  3306 

C  PH3  0375 

720  VT=0* 

11  =  1 
JN= J 


CHECK  FOR  ID  PROBLEM 


1 F  4  I MAX— 1)39,810,39 
SET  HOOP  STRESS  TO  ZERO 
810  DDVK-0. 

GO  TO  3306 

NOW  ME  MILL  SET  INDICES  FOR  A  HOOP 
STRESS  CALC.  OR  A  ECAL. 

39  IF1 J~ 119908,41, 47 

ME  ARE  IN  THE  BOTTOM  ROM 

41  KA=N 
KB=K 
V  V=l. 

UUU=1 • 

42  I F ( I— 1 i 9907 ,43,44 

ME  ARE  IN  THE  LOMER  LEFT  CORNER 

43  KL=K 
KR=K+1 
FD=1  • 

E=  1 . 

GO  TO  100 

44  I  FI  I— I MAX ) 46,45,9903 

ME  ARE  IN  THE  LOMER  RIGHT  CORNER 

45  KR=K 
KL=K-1 
E=l. 

FD=1. 

GO  TO  100 

ME  ARE  IN  BOTTOM  ROM,  BUT  NOT  AT 
AXIS  OR  RIGHT  BOUNDARY  OF  GRID. 

46  KR=K+1 
KL=K— 1 
E  =  1  • 


PH3  0385 
PH3  0390 
PH3  0395 
PH3  0400 
PH3  0405 
PH3  0410 
PH3  0425 
PH3  043GT 
PH3  0435 
PH3  0440 

PH3  0450 
PH3  0455 
PH3  0460 
PH3  0465 
PK3  0470 
PH3  0475 

PH3  0485 
PH3  0490 
PH3  0495 
PH3  0500 
PH3  0505 
PH3  0510 
PH3  0515 

PH3  0525 
PH3  0530 
PH3  0535 
PH3  0540 
PH3  0545 
PH3  0550 
PH3  0555 
PH3  0560. 
PH3  0565 


4 
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47 


46 


49 


50 


51 


52 


53 


54 


55 


59 


FD=1. 

PH3 

057C 

GO  TO  100 

FH3 

057 

PH3 

058 

NOT  IN  BOTTOM  ROW 

PH3 

056 

IF(J-JMAX)54,48,9903 

PH3 

05  9 1 

PH3 

059 

WE  ARE  IN  THE  TOP  ROW 

PH3 

060 

KA=K 

PH3 

060 1 

KB=K— IMAX 

UUU=1 o 

PH3 

06  It 

VV=lc 

PH3 

062( 

IF! I— 1)9907*50,51 

PH3 

062 

PH3 

0631 

WE  ARE  IN  UPPER  LEFT  CORNER 

PH3 

063 

KR=K+1 

PH3 

064C 

KL=K 

FD=1  • 

PH3 

064‘ 

E=l. 

PH3 

06  5 1 

GO  TO  ICO 

PH3 

066 

IF ( I—IMAX)53*52»9903 

PH3 

0661 

PH3 

067c 

WE  ARE  IN  UPPER  RIGHT  CORNER 

PH3 

067‘ 

KR=K 

PH3 

068  c 

KL-K— 1 

E  =  l. 

PH3 

068' 

FD=  1  - 

PH3 

069 1 

GO  TO  100 

PH3 

070( 

PH3 

070 

WE  ARE  AT  THE  TOP,  BUT  NOT  AT 

PH3 

07h 

THE  AXIS  OR  RIGHT  BOUNDARY  OF  GRID 

PH3 

071- 

KR=K+i 

PH  3 

07 2i 

KL=K~1 

PH3 

072 

E*l. 

PH3 

073 

FD=1. 

PH3 

073 

GO  TO  100 

PH3 

074i 

PH3 

074 

WE  ARE  NOT  AT  THE  TOP  OR  BOTTOM, 

PH3 

075 

CHECK  OUR  POSITION  WITH  RESPECT  TO  THE 

PH3 

075 

AXIS  ANO  RIGHT  BOUNDARY  OF  GRID. 

PH3 

076 

VV=lo 

PH3 

076 

KA=N 

PH3 

077 

KB=K— IMAX 

PH3 

077 

UUU=1. 

PH3 

078< 

I  FI  1-1)9907,59,56 

PH3 

078‘ 

PH3 

079 

WE  ARE  ALONG  THE  AXIS 

PH3 

079 

KL=K 

PH3 

080 

KR=K+1 

PH3 

080 

E=  1  o 

FD=  1  c 

PH3 

081 

on  non  ooono  on  o  on  no 
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GO  TO  100 

56  IFII-1MAX)58,57,9903 

WE  ARE  ALONG  THE  RIGHT  BOUNDARY  OF  CR10 

57  KR=K 
KL=K— 1 

E-i. 

FD=  1  • 

GO  TO  100 

WE  ARE  IN  THE  MESH,  NOT  AT  ANY  BOUNDARY 

58  E=i, 
f  D=  1  • 

KR=K+ 1 
KL=K— 1 

GO  TO  100 

100  IF 4  VT 1 102 , 101 ,102 

CALCULATE  THE  HOOP  STRESS  FOR  CELL  K. 

101  CALL  HOOP 
GO  TO  3306 

CALCULATE  THE  CHAN  2  IN  INTERNAL  ENERGY 
DUE  TO  THE  WORK  DONE  BY  THE  STRESSES. 

102  CALL  ECALC 
GO  TO  801 


CELL  K  IS  VOID,  SET  ALL  9 

STRESSES  TO  2ERO.  #****••**♦*♦***#**♦* 

3340  SNT=0. 

STT=0. 

SNL-0. 

SNfiaO. 

STL-0. 

STRsO. 

SNB=0. 

STB=0. 

D0VK-=0. 

XI =0. 

X2-0. 

GO  TO  3326 

RETURN  TO  HERE  AFTER  CALCULATING  THE 
HOOP  STRESS 

3306  1F( 1MAX— 1 19901,3311 *3310 

CHECK  FOR  KAREFIEO  MATERIAL  IN  THE  CELL  TO  THE 
RIGHT. 

3310  1F(AMX(K+11 ) 9902,3312,4000 

4000  IF ( AMX( K«-1J/(TAU4I*11*DY(J) )~AMDM*2( 1151)3312,3312,14 


PH3 

0820 

PH3 

0825 

PH3 

0830 

PH3 

0835 

PH3 

0840 

PH3 

0845 

PH3 

0855 

PH3 

0860 

PH3 

0865 

PH3 

0870 

PH3 

0875 

PH3 

0880 

PH3 

0885 

PH3 

0890 

PH3 

0895 

PH3 

0905 

PH3 

0910 

PH3 

0915 

PH3 

0920 

PH3 

0935 

PH3 

094Q 

PH3 

0945 

PH3 

0950 

PH3 

0955 

PH3 

0960 

PH3 

0965 

PH3 

0970 

PH3 

0975 

PH3 

0980 

PH3 

0985 

PH3 

0990 

PH3 

0995 

PH3 

1000 

PH3 

1005 

PH3 

1010 

PH3 

1015 

PH3 

1020 

PH3 

1025 

PH3 

1030 

PH3 

1035 

PH3  1045 


nror  noon  o  o  c>  o  o  o  on  non  non 
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3311 


3312 


14 


27 


32 


33 

723 

724 

725 


WE  ARE  AT  THE  RIGHT  BOUNDARY  OF  GRID 

PH3 

105 

SET  THE  STRESSES  ON  THE  RIGHT  =  TO 

PH3 

105 

THOSE  ON  THE  LEFT 

PH3 

106 

SNK=SNL 

PH3 

106 

STR=STL 

PH3 

107 

GO  TO  3316 

PH3 

107 

PH3 

106 

PH3 

108‘ 

THE  CELL  TO  THE  RIGHT  OF  CELL  K  IS  VOID. 

PH3 

109i 

SET  THF  STRESSES  ON  THE  RIGHT  TO  0. 

PH3 

1 09 

PH3 

1 10' 

SNR=0. 

PH3 

1101 

SIR=0. 

PH3 

lilt 

DDVK-0. 

S1=0. 

S2=0. 

S3=0. 

S4=0. 

S5=0. 

GO  TO  3316 

PH3 

ur 

PH3 

112C 

HE  ARE  PREPARING  TO  CALCULATE  THE  STRESSES 

PH3 

112" 

AT  THE  RIGHT  OF  THIS  CELL. 

PH3 

1 13c 

ME  ARE  NOT  AT  IMAX,  BUT  MAY  BE  AT  J=1 

PH3 

113" 

OR  J= JMAX  OR  IN  THE  MESH. 

PH3 

114 

PH3 

114" 

KR=K*1 

PH3 

115c 

11*1 

JN=J 

1F(J— 119908»2fl»27 

PH3 

115" 

PH3 

116( 

NOT  IN  BOTTOM  ROM 

PH3 

116 

IF< J-JMAXJ33, 32.9903 

PH3 

117c 

PH3 

117 

ME  ARE  IN  TCP  ROM 

PH3 

118t 

KA=K 

PH3 

118' 

KAR=KR 

PH3 

119( 

KB=K—IMAX 

PH3 

119 

KBR=KB+1 

PH3 

120t 

UUU=1. 

PH3 

120 

VV=1. 

GO  TO  31 

PH3 

1211 

PH3 

122 

ME  ARE  INSIDE  THE  MESH. 

PH3 

122 

NOTE.  THE  SPECIAL  BOUNDARY  CONDITIONS 
FOR  EMPTY  CELLS  ADJACENT  TO  CELL  K. 
KB=K- I MAX 
KBR=KB+ 1 

IF(AMX(KBI1 725.725# 726 
KB=K 
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• 

- 

726 

IF ( AMX4  KBRi i  727,727,730 

• 

727 

KBRsK+1 

73C 

VV=  1 . 

KA=N 

PH3 

1240 

KARSKA+1 

PH3 

1245 

GO  TO  31 

PH3 

1260 

C 

PH3 

1265 

C 

WE  ARE  IN  THE  BOTTOM  ROW. 

PH3 

1270 

28 

KB=K 

PH3 

1275 

C 

KBR-KR 

PH3 

1280 

c 

KA=N 

PH3 

1285 

c 

KAR=N+1 

PH3 

1290 

c 

UUUS1. 

c 

t/V=  1. 

PH3 

1300 

c 

31 

CALL  GRAOR 

PH3 

1305 

c 

C 

PH3 

1310 

c 

C 

CALCULATE  THE  VELOCITY  GRAOIENT S  AT  THE 

PH3 

1315 

c 

RIGHT 

PH3 

1320 

29 

r 

CALL  ST  RE SR 

PH3 

1325 

t 

c 

PH3 

1330 

c 

CALCULATE  THE  NORMAL  AND  SHEAR  STRESS 

PH3 

1335 

c 

c 

AT  THE  RIGHT. 

PH3 

1340 

*c 

30 

CONTINUE 

PH3 

1345 

c 

c 

PH3 

1350 

l 

c 

OONT  PUT  THE  INDIVIDUAL  STRESSES  INTO 

PH3 

1355 

« 

c 

ARRAYS  UNTIL  LATER. 

PH3 

1360 

GO  TO  3316 

PH3 

1365 

c 

PH3 

1370 

3316 

IF! JMAX-J 19903,3315,3320 

PH3 

1375 

C 

PH3 

1380 

C 

WE  ARE  AT  THE  TOP  OF  THE  GRID,  SET  THE 

PH3 

1385 

c 

STRESSES  AT  THE  TOP  *  THOSE  AT  BOTTOM  OF 

PH3 

1390 

c 

THE  CELL. 

PH3 

1395 

c 

PH3 

1400 

c 

PH3 

1405 

3313 

SNT-ASNU) 

PH3 

1410 

c 

STT*AST<I1 

PH3 

1415 

c 

D0VK*0. 

GO  TO  3325 

PH3 

1420 

C 

PH3 

1425 

C 

WE  ARE  NOT  AT  THE  TOP  OF  THE  GRID. 

PH3 

1430 

c 

PH3 

1440 

c 

CHECK  FOR  RAREFIED  MATERIAL  IN  CELL  ABOVE. 

3320 

IF1 AMXl N) 19906, 3322 ,4010 

6010 

IF(AMX(N)/(TAU(1 ) *DY l J+l ) )- \MDM*Z( 1 15 1) 3322,3322, 12 

( 

C 

PH3 

1445 

£ 

C 

THE  CELL  ABOVE  CELL  IK)  IS  VOID 

PH3 

1450 

c 

SET  THE  STRESSES  ABOVE  TO  0. 

PH3 

1455 

♦ 

3322 

SNT=0. 

PH3 

1460 

c 

» 

no  on  r o n o  nnnooooo 
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STT=0. 

ODVK=0. 

S6=G. 

S  7=0. 

S8=0. 

S9=0  . 
S10=0. 

GO  TO  3325 


kit  MILL  CALCULATE  THE  VELOCITY 
GRADIENTS  AND  STRESSES  AT  THE  TOP  OF 
THE  CELL. 

WE  ARE  NOT  AT  THE  TQP(J=JNAX1 

12  KA=N 
1 1  =  I 
JN=J 

21  IF ( I— il 990  7  « 15f 16 

16  IF(I~IMAXJ18,17,9901 

INSIDt  THE  MESH 

NOTE  THE  SPECIAL  BOUNOARY  CONDITIONS 
FOR  EMPTY  CELLS  ADJACENT  TO  CELL  K. 
18  KAR=KA+ 1 
KR=K*1 

732  IFUMXIKAR1 1733,733,734 

733  KAR-KA 

734  IF(AMXIKRJI735,735,736 

735  KR=K 

736  UUU=1. 

V  V= 1 . 

KL=K- 1 
KAL=K A- 1 
GO  TO  24 

ALONG  THE  RIGHT  BOUNDARY 

17  KAR=K A 
KR=K 

KAL=KA- 1 
KL-  K—  1 
UUU=1. 

VV=lo 
GO  TO  24 

WE  ARE  ALQM  THE  AXIS 
15  KL=K 
KAL=KA 


PH3  146 


PH3  147 
PH3  147b 
PH3  148( 
PH3  148 
PH3  149 
PH3  149 
PH3  150t 
PH3  150 
PH3  15  li 
PH3  151 


PH3  I52i 
PH3  152" 
PH3  153C 
PH3  153‘ 


PH3  154 
PH3  156 
PH3  156 
PH3  157 
PH3  157 
PH3  158 
PH3  158 
PH3  159 
PH3  159 
PH3  160 
PH3  160 

PH3  161 
PH3  162 
PH3  162 
PH3  163 
PH3  163 


o  o  o  n  o  o  n  o  n  non  nnnn  noon  o  no  non  on 


98 


CHECK  FOR  ID  PROBLEM 


IF( 1MAX-1 )811, 812,811 

812  KAR=KA 
KR-K 
UUU*1. 

GO  TO  813 
811  KAR=KA*1 
KRs'K+l 
UUU*i. 

813  VV=  1. 

CALCULATE  THE  VELOCITY  GRADIENTS  AT 
THE  TOP  OF  THE  CELL* 

24  CALL  GRADZ 

CALCULATE  THE  STRESSES. 

25  CALL  STRESZ 

26  CONTINUE 
GO  TO  6999 


ARRIVED  HERE  AFTER  COMPLETION  OF  THE 
CALCULATION  OF  THE  STRESSES  AT  THE  TOP. 
6999  IF (1-Jl 3325, 7001, 9908 


WE  ARE  IN  THE  BOTTOM  ROW,  NOW 
CHECK  THE  BOUNOARY  CONDITION  AT  THE  BOTTOM. 
7001  I~(CVI$ 17003, 7002, 7002 


BOTTOM  BOUNDARY  IS  REFLECTIVE 
7002  SNB*0. 

STB=0. 

GO  TO  3325 


BOTTOM  BOUNOARY  IS  TRANSMITT IVE,  SET  THE 
BOTTOM  STRESSES  TO  THE  TOP  (WE  JUST 
FINISHED  CALCULATING  THEM I • 

7003  SNB=SNT 
SIB=STT 
GO  TO  3325 

NOW,  WE  HAVE  ALL  THE  STRESSES  OF  CELL  K 
THUS  WE  CAN  CALCULATE  U  DOT  AND  VDOT. 
3325  CONTINUE 


( 

< 


PH3  1645  ( 

( 

PH3  1660 
PH3  1665 
PH3  1670 
PH 3  1675 
PH3  1680 
PH3  1685  ( 

PH3  1690  ( 

( 

PH3  1695  ( 

PH3  1700 
PH3  17C5 
PH3  1710  *< 

PH3  1715  ( 

PH3  1720  ( 

PH3  172*5  . < 

PH3  1730  ( 

PH3  1735  ( 

PH3  1740  ( 

PH3  1745 
PH3  1750 
PH3  1755 
PH3  1760 
PH3  1765 
PH3  1770  ( 

PH3  1775  ( 

PH3  1780  ( 

PH3  1785  ( 

PH3  1790  ( 

PH3  1795 
PH3  1800 
PH3  1805 
PH3  1810  ( 

PH3  1815  ( 

PH3  1820 
PH3  1825  < 

PH3  1830  ( 

PH3  1835  ( 

PH3  1840 
PH3  1845  ( 


non  o  o  ooooo  r>  r>  r>  r>  n  rn  noon  n  o  no 
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300  CONTINUE  PH3  lb 

JN=  J 
1 1  =  I 

PH3  lb 

CALL  DELTAU! COMPUTES  UDQT)  PH3  16 

IF! J-l)310,310,311  PH3  lb 

311  STB=AST ( I J  PH 3  16 

310  CALL  DELTAU  PH3  18 

PH3  18 

30.1  CONTINUE  PH3  18 

NOW  CALL  DELTAV! COMPUTES  VDOT)  PH3  18 

302  CONTINUE  PH3  .18 

I F 1 J —  1 i 312,312,313  PH3  19 

313  SNB=ASN( I )  PH3  19 

312  CALL  DELTAV  PH3  19 

303  CONTINUE  PH3  19 

PH3  19 

BY  NOW  WE  HAVE  THE  ACCELERATION  PH3  » 

(BOTH  COMPONENTS)  OF  CELL  K  DUE  TO  PH3  lv 

THE  STRESSES*  PH3  19 

304  IF(i-l)9907, 3326, 305  PH3  19 

305  I F( I— I MAX) 306,3326,9901  PH3  19 

PH3  19 

/////////////✓//////✓/////////////////////////// 

CHECK  FOR  OVERSHOOT  IN  THE  RADIAL  DIRECTION.  PH3  19 


AT  THE  LEFT  INTERFACE  OF  CELLIKJ, 
////////////✓//////////////////////////////////✓ 
CALCULATE  DELTA  U  AT  CYCLE  N/ 

DELTA  UDOT  AT  CYCLE  N+l 

306  WS=— 1  Ul  K)— U ( K-l J )/ I XI— DUDOT ! I— 1 1 ) 

307  IFIWSJ450, 308,308 

308  IF( DT—W SI 450,309, 309 
450  UPR=0. 

GO  TO  400 

CHECK  IF  WS  IS  BETWEEN  0.  AND  DT 
IF  LESS  THAN  ZERO,  BYPASS  CHECK. 

IF  IT  IS  LESS  THAN  DT,  REDUCE  THE 
STRESS! NORMAL) BY  THE  RATIO  WS/DT. 

309  WS=WS/DT 
SNL=SNL*WS 
UPR-—  1® 

CHECK  THE  OTHER  COMPONENT. 

400  WS=— (  V!K)— V(K— 1))  /  (  X2— DVDQT  (  I-D) 

NOW  CHECK  THE  SIGN  AND  MAGNITUDE  WITH 
RESPECT  TO  DTIHYDRO). 

*  401  IF(WS) 501,408,408 

C 


PH3  19 
PH3  19 
PH3  19 
PH3  19 
PH3  19 
PH3  19 
PH3  19 
PH3  19 
PH3  20 
PH3  20 
PH3  20 
PH3  20 
PH3  20 
PH3  20 
PH3  20 
PH3  20 
PH3  20 
PH3  20 
PH3  20 
PH3  20 
PH3  20 
PH3  2C 
PH3  20 


% 


non  n  o  nooooooon  o  o  o  on  on 


100 


IFI-1 J8YPASS  CHECK 

IF  GREATER  THAN  0.  CHECK  AGAINST  0T 

408  I FI  DF— HSJ  501*409*409 

IF! GREATER! BYPASS  CHECK 
501  UPZ=0. 

GO  TO  500 

409  WS-HS/DT 

REDUCE  SHEAR  STRESS  AT  THE  RIGHT 

BY  WS/DT 

STL=STL*WS 

UPZ=-1. 

IF  PIK1=0.  NO  OVERSHOOT  IN  EITHER 
COMPONENT. 

IF  P I  K )  =—  1  •  THE  SHEAR  STRESS  HAS 
MODIFIED. 

IF  P I Ki =1*  THE  NORMAL  STRESS  HAS 
MODIFIED. 

IF  PI K)=2.  BOTH  OF  THE  STRESSES 
REQUIRED  MODIFICATIONS. 

500  CONTINUE 

I F I UPR1 960 1 , 9602 , 9602 

9601  IFI UPZ 19603 *9604 *9604 

9603  PIK—1 J=2. 

GO  TO  3326 

9602  IFI UPZJ 9605*  3326*3326 
9605  P I K— 1 J  =— 1 . 

GO  TO  3326 

9604  PIK—1 1=1.0 
GO  TO  3326 

******  NOTE*  WE  ONLY  SET  THE  COMPLETE  ARRAY  FOR  THE 
FIRST  ROW  ************ 

3326  IFI J-l) 9902*3328*601 
3328  ASNI I 1=SNF 
ASTI  I )  =  STT 
RSNI I ) =SNL 

640  RSNI I+I )=SNR 
RST 1 1*1 )=STR 

641  ASNBI I ) =SNB 
RST I Ii=STL 
ASTBI I )=STB 
SIG33I I )=DDVK 
DUDOT I i 1=X1 
DVDOT ( 1 )=X2 


PH3  2080. 
PH3  2085 
PH3  2090 
PH3  2095 
PH3  2100 
PH3  2105 
PH3  2110 
PH3  2115 
PH3  2120 
PH3  2125 
PH3  2130 
PH3  2135 
PH3  2140 
PH3  2145 


Ph3  2435 
PH3  2440 
PH3  2445 
PH3  2450 
PH3  2455 
PH3  2460 
PH3  2465 
PH3  2470 
PH3  2475 
PH3  2480 
PH3  2485 
PH3  2490 
PH3  2495 
PH3  2500 


SET  THE  RIGHT  STRESSES  FROM  CELL  IKJ  TO 
=  THE  LEFT  STRESSES  FOR  CELL  IK+1J. 
SNL=SNR 


PH3  2505 


anno  o  n  n  o  ana  no  non  non  onoooooono 
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STL-STR  PH3  25 

GO  TO  3360 

PH3  25 

//////////////////////////////////////////////// 

CHECK  FOR  OVERSHOOT  DUE  TO  STRESSES  PH3  25 

AT  THE  TOP  OF  CELL ( K- IMAX I  OR 
THE  BOTTOM  OF  CELL ( Ki  • 

////////////////////////////////////////////✓//// 

*********  PH3  25 

CAME  HERE  AFTER  CALCULATING  A  CELL  IN  A  PH3  25 

ROW  DIFFERENT  THAN  FOR  J=1  .  PH3  25 

**********  PH3  25 

601  KRM=K-I MAX  PH3  25 

IF( AMX( KRM ) / I TAU1 1 1 *DY( J-i) )-AMDM*ZU 15 J ) 901 , 901, 602 

**********  PH3  25 

KRM= INDEX  OF  CELL  BELOW  '*  PH 3  25 

**********  PH3  25 

602  WS=- ( UI  K|— U< KRM) 1/  ( Xl—DUDQT CUT  PH3  25 

**********  PH3  25 

REMEMBER,  THE  STRESSES, UOOT, VDOT,  ANO  PH3  25 

HOOP  STRESS  FOR  CELL  K  HAVE  NOT  BEEN  PLACED  PH3  25 

IN  THE  ARRAY.  PH3  25 

*******4 ***  PH3  26 

603  IF(WS)604,605,605  PH3  26 

605  IFIDT-WS)604,606,606  PH3  26 

604  UPR=0.  PH3  26 

GO  TO  607  PH3  26 

**********  PH3  26 

REDUCE  THE  SHEAR  STRESS  ON  THE  TOP  PH3  26 

**********  PH3  to 

606  WS— WS/DT  °H3  26 

ASTI  IJ  =  ASTU)*WS  FH3  26 

UPR=— 1.  PH3  26 

CHECK  OTHER  COMPONENT  PH3  26 

607  WS=-I  VlK»-VIKRMn/(X2-OVDOTI  I)  )  PH3  26 

608  IF(WS)61«, 615,615  PH3  26 

615  IF(DT—W 51614,616,616  PH3  26 

614  UPZ=Q.  PH3  26 

GO  TO  617  PH3  26 

***********  PH3  26 

REDUCE  THE  NORMAL  STRESS  ON  THE  TOP  PH3  26 

**********  PH3  26 

616  WS=WS/DT  PH3  27 

ASNII)-ASN(I) *WS  PH3  27 

UPZ=-1.  PH3  27 

***********  PH3  27 

SET  INDICES  FOR  CELL  BELOW  FOR  POSSIBLE  UPDATING  PH3  27 

OF  UDOT  AND  VDOT  PH3  27 

***********  PH3  27 

617  KK-K  '  PH3  27 


o  o  o 
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KKK=N 
K-K-“  I  "'AX 
N=KK 
11=1 
JN=J " 1 

***$$«*$*****$* 

SAVE  THE  STRESSES  OF  CELL  K 
TEMPORARY. 

SNB=ASN1 1 ) 

STB=ASTU) 

PK  ( 1 1 =SNB 
PM 2)=S1B 
PM3)  =  STR 
PM4J  =  SNR 
PK ( 5 ) =SNT 
PK16)=STT 
PK ( 7 ) =SNL 
PM 8J=STL 
PK( 9 j=ODVK 
PK ( 10J=X1 
PK i 1 1 )  =  X2 

NOW  SET  STRESSES  FROM  K-IMAX  INTO 
THE  CORRECT  STORAGE  FOP.  THE 
SUB— ROUTINES. 

626  SN8=ASNB( I ) 

STB=ASTBUJ 
STR=RST<I+1J 
SNk=RSNU  +  l) 

SNT=ASN ( I J 
STT=ASTII) 

SNL=RSN4 I ) 

STL=RSTU) 

DDVK=SI G33 (  I  i 
X1=DUDGT1 Ii 
X2=DVDOT(I) 

620  IFUJP2) 624,9700,9700 

624  CALL  DELTAV  - 

GO  TO  9702 

9700  IFlP(K)-lo )9710, 9701.9710 
9710  IFIP(K) 1624,9701*624 

9701  CONTINUE 

9702  IF { UPR 19703*9704 *9704 

9703  CALL  DELTAU 
GO  TO  9705 

9704  IF(P(K)~1. >9705,9703, 9703 

9705  CONTINUE 

625  IF(J— 2) 800,650,651 
650  SIGC ( 1 1 =U<  K ) 

GAMcm=vm 


PH3  2740. 
PH3  2745 
PH3  2750 
PH3  2755 
PH3  2760 
PH3  2765 
PH3  2770 
PH3  2775 
PH3  2780 
PH3  2785 
PH3  2790 
PH3  2795 
PH3  2800 
PH3  2805 
PH3  2810 
PH3  2815 
PH3  2820 
PH3  2825 
PH3  2830 
PH3  2835 
PH3  2840 
PH3  2845 
PH3  2850 
PH3  2855 
PH3  2860 
PH3  2865 
PH3  2870 
PH3  2-75 
PH3  2880 
PH3  2885 
PH3  2890 
PH3  2895 
PH3  2900 
PH3  2905 
PH3  2910 
PH3  2915 


on on  n  nno  on  on 
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GO  TO  800 

651  KBB=K-I  MAX 
WSU=UI KBB) 

WSV=V ( KBB ) 

UlKSBJ  =  SIGCUJ 
V(KBB)=GaMC<  IJ 

SET  FLAG  TO  1.IUSE  THE  LOOK-UP 
ROUTINE  THAT  THE  HOOP  STRESS  USES) 

800  IF( i— 2 ) 802*  603*803 

SET  THE  VELOCITIES  AT  THE  LEFT  AND  SAVE 
THE  OLD  ONES. 

803  UTEF=UI K-l ) 

VTEF=VlK-l) 

UI K— 1 )  =  Si GC ( 1-1 ) 

VI K—l ) =GAMC 1 1-1 ) 

802  VT= 1« 

GO  TO  39 

**  FINALLY  GETTING  AROUNO  TO  CALCULATE 
THE  VELOCITIES  AND  INTERNAL  ENERGIES. 
CALCULATE  THE  RADIAL  COMPONENT,  AND  AXIAL 

801  AIXIK)  =  AIX4KH-VISC 
670  CONTINUE 

IFIAIXI KJ—VVABOV) 700,700,656 
700  SUM=SUM*AIXIKI*AMX«K) 

A IX  I K ) =0. 

656  IF! 1—2)652,805,805 

RES?1’  THEIN  +  1)  VELOCITIES  IN  CELL  IK-1) 
805  U ( K— 1 ) -UTEF 
V(K-1)=VTEF 
GO  TO  652 

652  IFIJ— 3)653,654,654 
654  SIGCm=UIK) 

GAMCI I )=V(K) 

UIKBBi=WSU 
VIKBB)=WSV 
GO  TO  653 


PH3  2995 
PH3  3000 


PH3  3010 
PH3  3015 
PH3  3020 
PH3  3025 


/////////////////////////////////////// ///////✓//✓///////// 
NOTE,  HERE  WE  CHECK  ON  THE 
POSSIBLE  OVERSHOOT  FROM  THE  HOOP  STRESS 
653  WSA=X1*DT 
WS=UI K) +WSA 
IFIU(K) 1661,658,657 

657  IF(WS)660,658,658 

658  UIKl-WS 
GO  TO  659 

660  DX1=2.*PI0Y/AMX(K)*DX1I)*DY( JNi 
DX1=— 0X1*DDVK 
X1=X1-0X1 

LXi  =  DXl*ABS(lj(K)/WSA) 


onoriooooooonoooo  r  no 
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661 

659 

701 

3 

3  702 

703 


70s 


873 


900 


901 


X1=X1+DX1 
U( KJ=U(  Kl+Xi*DT 
GO  TO  659 
IF(WSJ658, 658,660 
V(K)  =  V<KJ<-K2*0T 

I F  ( A8  S  t  U  ( K ) 1-VVBLQ1 701»  701 1  702 

SUM=SUM«-U(K1**2/2.*AMX(K» 

U(K)=0. 

I F ( ABSi V( K ) )— VVBL0) 703,703,704 
SUM=SUM*V  ( K )  **2/2.*AMX( K) 
ViK)=0. 

CONTINUE 


Vise  IS  DI/OT 

PH3 

3090 

***  RESET  THE  INDICES 

NOW 

,  K 

AND  N 

PH3 

3095 

WERE  TEMPORARY  SET  FOR  CELL  BELOW 

PH3 

3100 

K=KK 

N=KKK 

PH3 

3110 

***** ***** ** 

PH3 

3115 

PH3 

3120 

WE  NOW  HAVE  COMPLETED 

THE 

INTEGRATION 

OF 

PH3 

3125 

THE  MOMENTUM  AND 

ENERGY  EQUATIONS 

PH3 

3130 

FOR  CELL  K-IMAX 

PH3 

3135 

PH3 

314  6 

*********** 

PH3 

3145 

WE  CAN  NOT  PLACE 

ALL 

THE 

STRESSES 

OF 

CELL 

K 

PH3 

3150 

INTO  THE  (11 

PH3 

315S 

ARRAY,  SINCE  THE 

OVERSHOOTiON 

THE 

TOPI 

OF 

CELL 

PH3 

3160 

K-IMAX+1  HAS  NOT 

BEEN 

CHECKED 

PH3 

3165 

PH3 

3170 

SET  THE  STRESSES 

THAT 

WE 

STORED 

PH3 

3175 

TEMPORARY  IN  PK(l)  THRU  PK(9J 

,  BACK 

INTO 

PH3 

3180 

THE  STORAGE  THAT 

THE 

SUBROUTINES 

PH3 

3185 

RECOGNIZE. 

PH3 

3190 

SNB=PKI 11 

PH3 

3195 

STB=PK( 2) 

PH3 

3200 

STR=PK(3J 

PH3 

3205 

SNR=PK<4) 

PH3 

3210 

SNT=PK( 5) 

PH  3 

3215 

STT=PK(6J 

PH3 

3220 

SNL=PK( 7 ) 

PH3 

3225 

STL=PK( 8) 

PH3 

3230 

DDVK=PK ( 9 ) 

PH3 

3235 

Xi=PK( 10) 

PH3 

3240 

X2=PK(111 

PH3 

3245 

ASNB ( I 1 =ASN ( I i 

PH3 

3250 

ASTB( I ) =AST ( I ) 

PH3 

3255 

RSN( I )=SNL 

PH3 

3260 

RST  (  I  )-=STL 

PH3 

3265 

ASN( I 1 =SNT 

PH3 

3270 

ASK  I  )  =  STT 

PH3 

3275 

UipUIOUtOUlOUlUWUVIUUUJUiUWiwwnju'wifwv1  w»v> 
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.c. 


650 

851 

860 

661 

C 

C 

C 

902 


>  904 

► 

>  C 

998 

) 

>  999 

>  C 

>  C 

)  903 

• 

599 


C 

C 

C 

3360 

3361 
C 

C 


952 

953 
C 

3302 

C 

c 

c 

c 


906 

907 


UPDATE  CELL  K  IF  STRESS  WAS  CHANGED 
If (UPZ)850,b51,851 
CALL  DELTAV 
IF ( UPR) 860 ,861,861 
CALL  DELTAU 
CONTINUE 
♦♦♦♦♦♦♦*♦*♦ 

NOW  SET  THE  UDOT  AND  VDOT  AND 
HOOP  STRESS  INTO  THE  PROPER  ARRAY. 

DUDOT (I )=X1 
DVDOT ( 1  |=X2 
S I G33 ( I J=DDVK 
CONTINUE 

I F  4  1-2)999,998,998 

HERE,  SET  THE  STRESS  ON  THE  RIGHT  INTO  THE  (II  ARRAY 
RSNI I )=SNL 
RST ( I )=STL 
CONTINUE 

SET  THE  STRESS  AT  THE  LEFT  OF  CELL  <K*1>  TO  THE 

RIGHT  OF  CELL  4  Ki 

SNL=SNR 

STL=STR 

I F ( I— IMAX )  3360,599,599 
RSNU*1)  =  SNR 
RST ( 1+1)  =  STR 
GO  TO  3360 

*****  END  OF  DO  LOOP  ON  (I)  ****** 

K=K+1 
N=N+ 1 
CONTINUE 

CHECK  HERE  AT  THE  RIGHT  OF  ACTIVE 
GRID. 

KLAST=K— 1 

IF  lABS(UIKLAST) )+ABS (V(KLAST)J +A IX (KLAST ) )  952,953,952 

NRC=1 

CONTINUE 

*****  END  OF  DO  LOOP  IN  THE  4 Ji  DIRECTION  ************ 
CONTINUE 

HERE  WE  WILL  UPDATE  THE  LAST 
ROW. 

SAVE  THE  OLD  VELOCITIES 

HERE  WE  NEED  THE  OLD  VELOCITIES 

K  =  ( 12—1 )*IMAX+2 

JN=I  2 

N=K— IMAX 

DO  980  1=1,11 

IF( I— 1)908,907, 908 

KL=K 


PH3  328, 
PH3  328b 
PH3  329C 
PH3  329b 
PH3  3301 
PH3  330b 
PH3  3330 
PH3  333b 
PH3  3340 
PH3  334b 
PH3  3350 
PH3  335b 
PH3  3310 
PH3  331b 
PH3  332 C 
PH3  332b 


PH3  37lb 
PH3  3720 
PH3  372b 
PH3  3730 
Ph3  373b 
PH3  3740 


PH3  374b 
PH3  3750 


o  r>  n  o  r»  r>  on 
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GO  TO  909 

908  KL=K“ 1 

909  IF 1 I2-JMAX)91i,910,911 

910  KA=K 
GO  TO  912 

911  KA=K*IMAX 

912  IFII-IM  AXJ914, 913,914 

913  KR=K 
GO  TO  915 

914  KR=K+1 

NOW  WE  HAVE  THE  INDICES  FOR 
SUBROUTINE  ECALC* 

915  IF! I— 11917*916*917 

916  FLEFT ( 1 J=U(K ) 

YAMC ( I J  —  V €  K J 
WSU=U(N1 
WSV=V(Nj» 

VELOCITIES! PHASE  1)  FROM  CELLS  BELOW 
ARE  AVAILABLE  FROM  THE  MAIN  LOOP 
AND  ARE  STORED  IN  THE  ARRAYS  1 S IGC  AND  GAMC  ). 
U!N)=SIGC(IJ 
V! N)=GAMC( I i 
GO  TO  919 

917  WSU=U(N) 

WSV=V4N) 

UTEF=U(K— 11 
VTEF=V!K-11 

SAVE  THE  UPDATED  VELOCITIES 

918  U(NJ  =  S1GCUJ 
VIN)*GAMC( I ) 

U!K— 1)=FLEFT(I— 11 
V!K— 1)=YAMCI I— 1) 

919  II  =  I 
KB-=N 

SET  THE  STRESSES  FROM  THE  ARRAYS 
INTO  THE  SINGLE  STORAGE. 

SNT=ASN1 I ) 

STT=AST ( I ) 

STR=RST !  I  +  ll 
SNR=RSN( 1+1) 

SNB= ASNBC I ) 

STB=ASTBCIJ 
S..  =RSN(  I J 
STL=RST ( 1 1 
Xl=OUDOT( II 
X2=DVDOT(U 
ODVK=SIG33(  1 1 
IF ( JN— UMAX ) 530,940  *930 
C  IF  WE  ARE  AT  THE  TOP  BOUNDARY  OF 

C  THE  GRID, SET  THE  STRESS  GRADIENTS 


on  o  o  o 
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C  TO  ZERO 

940  SNI-3NB 
STT=STB 

930  IF( I“IMAX)943,941,943 

If  RE  ARE  AT  THE  RIGHT  BOUNDARY 
OF  THE  GRID,  SET  THE  STRSS  GRADIENTS 
TO  ZERO 

941  STR=STL 
SNR=SNL 

943  CALL  DfcLTAU 
CALL  DELTAV 

920  CALL  ECALC 

921  AIX(KJ=AIX(K!+VISC 
991  CONTINUE 
994  CONTINUE 

IF ( Ai X( K )- VVABOV 1705,705,922 

705  SUM=SUM+AIX(KJ^AMX1K) 

A IX ( K  )  =  0. 

922  FLEFT ( I 1=U( K ) 

YAMC ( I ) =V( K) 

923  U  (  K )  -U(  K  J  #-DT*Xl 
V1K)=VIKJ*DT*X2 

IF ( ABS( 0( K ) )— VVBLO) 706,706,707 

706  SUM=S0M*U(K)**2/2.*AMXlKl 
U(K)=0. 

707  IF I AB SI V(KJ)-VVBLO) 708,708,709 

708  SUM=SUM+V(K)**2/2.*AMX|K1 
V(K)=0. 

709  CONTINUE 
951  IF(I-l)  925,926,925 
926  U(Ni=WSU 

V ( N ) =WS V 
GO  TO  924 

RESET  THE  NEW  VELOCITIES  FOR  THE  LEFT 
AND  BOTTOM  CELLS. 

925  U ( N)  =  WSU 
V ( N )  =  WSV 
U ( K— 1 J =UTEF 
V ( K-l ) =VTEF 
GO  TO  924 

924  K=K+ 1 
N=N«-i 

980  CONTINUE 
C  CHECK  HERE  AT  THE  TOP  OF  ACTIVE 
C  MESH. 

K  =  K~"  1 

IF(ABS(U1KJ1*ABS(VIK))*AIX(K)J 950,954,950 
950  NRT=1 
954  CONTINUE 

C  NOW  INCREASE  ACTIVE  GRID  COUNTERS  IF 
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C  NEEDED. 

21  =  2 1+NRC 
12=  I2+NRT 

IF(  2 1—1MAX 16100,6100, 6200 


6200 

2 1=  IMAX 

62  00 

IF ( I 2- J MAX 16201,6201*6202 

6202 

12=JMAX 

6201 

CONTINUE 

GO  TO  7777 

PH3 

3765 

9908 

NK=39 

PH3 

3770 

GO  TO  9999 

PH3 

3775 

9907 

NK=42 

PH3 

3780 

GO  TO  9999 

PH3 

3785 

9903 

NK=44 

PH3 

3790 

GO  TO  9999 

PH3 

3795 

9901 

NK=3306 

PH3 

3800 

GO  TO  9999 

PH3 

3805 

9902 

NK=3310 

PH3 

3810 

GO  TO  9999 

PH3 

3815 

9904 

NK=3320 

PH3 

3820 

GO  TO  9999 

PH3 

3825 

9900 

NK=3305 

PH3 

3830 

GO  TO  9999 

PH3 

3835 

9999 

NR=123 

PH3 

3840 

CALL  DUMP 

PH3 

3845 

7777 

SUMX=0. 

DO  9001  1=1,21 

K=  I  + 1 

00  9002  J=1 , 12 

P  ( K 1 =0. 

IF (A  2X1 K) )71C, 9002, 9002 

» 

710 

SUMX=SUMX«-AIXIK)*AMX(K1 

AIX(K1=0. 

9002 

K=K*IMAX 

9001 

CONTINUE 

ETH=E  TH-SUH— SUHX 

9000 

CONTINUE 

VT=VSAVE 

DT=DTI 

7778 

RETURN 

END  PH3  3855 
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•WBFTC  GRADR  LIST, DECK, REF 
SUBROUTINE  GRADR 

C  CALCULATES  THE  VELOCITY  GRADIENTS  ON 
C  THE  RIGHT  SIDE  OF  THE  CELL 

C  GRADR  REQUIRES  THE  FOLLOWING  1NOICES. 

C  I  , J,K,KR,KA ,KAR,KB,KBR« 

C  S1=DU/DR*** 

S1=(UIKR)-UIK1  l/DXiin 
1FIABSI  S1T-Z(107m,l,2 

1  S 1=0. 

2  CONTINUE 

C  S2=DV/DR*** 

S2=IV(KRJ— VIK))/DXII I) 

I F ( ABSI  S2)— 2 ( 107 ) ) 3, 3, A 

3  S2=0« 

4  CONTINUE 
WS=2-*DY( JN) 

C  S3=DU/DZ*** 

S3= I ( U( KA) +UIKAR) )/2«-(UlK8)+U( KBR1 )/2« I/WS 
IF<ABS<S3)-ZI107) J5,5,6 

5  S3=0o 

C  6  CONTINUE 

C  S4-0V/DZ*** 

S4=  I { VI KA ) +V ( KAR) ) /2«  *VV— I  V I  KB ) ♦VI KBR) ) /2. *UUU) /WS 
IFIABSI SAl-ZI 107)) 7,7,8 
*  7  S4*0- 

6  CONTINUE 

C  S5=U/R*** 

IFIGAMJS, 12,9 

12  S5=(UlKR)+U(Ki)/I2«*XiiI) ) 

IF(ABSIS5)-Z(107))9,9,10 
9  S5=0» 

10  CONTINUE 
RETURN 
END 


$ IBFTC  GRADZ  LIST, DECK, REF 

SUBROUTINE  GRADZ 

CALCULATES  THE  VELOCITY  GRADIENTS  AT  THE 
TOP  OF  THE  CELL 
S6=i)U/DZ*** 

S6=IUIKA)-UIK) )/DY( JN1 
IFIABSI S61-ZI 107) ) 1,1,2 
1  S6=0o 
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2  CONTINUE 

C  S7-DV/D2*** 

S7-( V(KA)-V{KJ )/OY ( JN) 

IF(ABS(S7)-Z(lC7J)3t3»4 

3  S7-»0« 

4  CONTINUE 
NS=2.  *OX<  in 

C  S8=DU/DR*** 

S8=( ( U< KARJ +UCKR) )/2.*VV-iU( KAL)*U(KL)1/2.*UUU1 /WS 
IF(ABS(  S8)-2(107))5,5.6 

5  S8=0. 

6  CONTINUE 

C  S9=D V/DR  ************ 

S9=I I  V  < KARJ *VIKR))/2«-(V{KALl*V{KL))/2« )/WS 
IFI ABS( S9J— 2 (107i)7»7»8 

7  S9-0. 

8  CONTINUE 

C  S 10=U/R*** 

IF<GAH)9,12.9 

12  S10=(U(KA)^U(Ki)/(X(II)4-XUI-m 
IF( ABS1 S10)-Z( 107) )9,9, 10 

9  S10=0. 

10  CONTINUE 

RETURN 

END 
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.  *IBFTC  STRESR  LIST  »DECK ,REF 
SUBROUTINE  STRESR 

C  THIS  ROUTINE  CALCULATES  THE  NORMAL 
C  AND  SHEAR  STRESS  AT  THE  RIGHT 
HAND  BOUNDARY  OF  THE  CELL 
CALCULATE  SNR4 SIGMA  II  AT  THE  RIGHT) 

NOTE,  ME  CAN  HAVE  VISCOSITY  OR  STRENGTH  OR  BOTH 

NEWT=U 

IF (DKE) 1,1,91 
91  If (DDXN)93,93,1 
93  NEWT-1 
B=0. 

GO  TO  100 

1  WS=.66666*ISl*SU-S4*S4+S5*S5)+.5* 

1 I S3+S2 ) **2 

THAT  MAS  THE  STRESS  DEVI ATQR= 

EPSILON  DOT  IAB)  *  EPSILON  OOT  <AB) 

K  IS  STORED  IN  DDXN 
MSA=WS*DX( ili*DX(ll) 

IF1SQRT (MSA)— Z( 112)*DXN)3,3, 4 

3  B=0. 

GO  TO  100 

4  CONTINUE 
MSAsBlG 

NOTE,  BIG  IS  CALCULATED  IN  PH3 
IF1 ABSI SI)— MSA) 10,10,11 

11  BUGR=ODXN 
GO  TO  12 

10  BUGR-ODXN*ABSI SI ) /MSA 

12  B-SQRT (2«*BUGR*BUGR/WS) 

100  ED0T1 l-Sl 

NOTE,  S1=DU/DR 
S1=0U/DR,  S4=OV/OZ,  S5=U/R 
EDQTAA=S1«-S4*S5 
EPDi 1=EDQT11— EDOTAA/3* 

NOW  CALCULATE  THE  NORMAL  STRESS 
SNR= ( B+DKE) *EPD11 
NOW  THE  SHEAR  STRESS 
DELTA  12  IS  ZERO,  THUS 
EPD12=ED0T12 
ED0T12=IS3+S2)/2. 

C  S3=DU/DZ,  S2=DV/DR 

IF(NEWT)95, 95,96 
96  B=0. 

GO  TO  97 

95  IF (ABSI S4)— WSA) 13,13,14 
14  BUGZ=DDXN 
GO  TO  15 

13  BUGZ=DDXN*ABSI S4)/WSA 
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15  B=SQRI(2.*BUGZ*BUGZ/WS) 
97  STR=( B+DKE ) *£D0T12 
RETURN 
ENO 


$ I8FTC  STRESZ  LIST .DECK, REF 
SUBROUTINE  STRESZ 

THIS  ROUTINE  CALCULATES  THE  NORMAL  AND  SHEAR 
STRESS  AT  THE  TOP  OF  THE  CELL. 

CALCULATE  SNT (SIGMA  22  AT  THE  TOP) 

NOTE,  WE  CAN  HAVE  VISCOSITY  OR  STRENGTH  OR  BOTH 
NEWT-0 

IF (DKE) 1,1.91 
91  IF (DDXNJ93.93 ,1 
93  NEWT=1 
B=0. 

GO  TO  100 

1  WS=«66666*(S8*S8+S7*S7«-S10*S10H-.5* 

1(S6+S9)**2 
WSA=WS*DY( JNJ*DYI JN) 

I F I  SORT ( WSA1— 21 112)*0XN)3,3,4 

3  B=0. 

GO  TO  100 

4  CONTINUE 

C  NOTE  . S7=DV/DZ 

WSA=BIG 

ir-IABSIS7)-NSA)6,6,7 

7  BUGR=ODXN 
GO  TO  8 

6  BUGR=0DXN*ABSIS7)/WSA 

8  CONTINUE 

B=SQRT(2.*BUGR*BUGR/WS) 

C  NOW  SIGMA  22  =  B*EP022 

100  £PD22=S7— ( S8+S7+S10) /3. 

SNT=| B+OKE I *EPD22 
EPD21=(S6+S9)Z2. 

C  NOTE,  S8=DV/DR 

IF! NEWT 1 95 , 95 ,96 
96  B-0. 

GO  TO  97 

95  IFIABSI S8J-WSA) 10,10,11 
11  BUGZ=DDXN 
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GC  TO  12 

10  BUG2=DDXN*ABS( S8J/WSA 
12  B=SURT 1 2«*BUGZ*BUGZ/RS) 
97  STT=(8+DKE)*EPD21 
RETURN 
END 
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$ 13FTC  HOOP  LIST, DECK, REF 

SUBROUTINE  HOOP 

HERE  WE  WILL  CALCULATE  THE  HOOP 
STRESS  FOR  CELL  K. 

K  IS  IN  (DDXNJ 

ETA  ZERGIVISCOSITY1IS  IN  DKE 
THE  HOOP  STRESS  IS  STORED  IN  CDVK 
GAM  IS  A  FLAG  FOR  THE  TYPE  OF 
COORDINATE  SYSTEM. 

IF (GAM) 101,102,101 

101  DDVK=0. 

GO  TO  103 

102  WS=XUI  J  +  X(  11-11 
EDQT33=UIKJ/WS*2. 

IF (DDXNJ 1,1,2 

C  WE  ASSUME  WF  HAVE  A  VISCOUS  MATERIAL. 

1  B=DKE 

GO  TO  100 

C  CALCULATE  8, WE  HAVE  A  RIGID  PLASTIC  MATERIAL 

2  WSR=DX(  III 
WSZ=DY  J  JN1 

WSA=( <UIKR1*E-U(KLJ*FD)/(2.*WSR1 1**2 
IF! ABSIWSAl-Z (107) 110,10,11 

10  WSA=0 • 

11  CONTINUE 

WSB=( (V(KA1*UUU-VV*V(KB11/I2**WS21 1**2 
IF (A8S(WSBJ-Z( 107)112,12, 13 

12  WSB=0. 

13  CONTINUE 
WSC=(2.*U(K1/WSJ**2 

WSD=( U<  :iAJ— jIKBI  )/(2.*WSZl 

'•  SD=(  (V4KR1-V#KL)  1/(2. *WSRUWSD  1**2 

IFIABSt  WSD1— ZI107) 114,14,15 

14  WSD=0. 

1C  CONTINUE 

B=. 60666* ( WSA*WSB*WSC1+.5*WSD 
WSA= ( OX ( 1 1 ) +DY ( JN) 1 /2  * 

WSA~B*WSA*WSA 

IF1SQRT1 WSA1— Z(1121*DXN13,3,4 

3  B=0. 

GO  TO  100 

4  CONTINUE 

B=SQRT(2.*DCXN*DDXN/8) 

C  NOW  SIGMA  33-B*EDQT  OF  33 

100  D0VK=B*£D0T33 

103  RETURN 
END 
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<  4 iBFTC  DELTAU  L I  ST .DECK , REF 
SUBROUTINE  DELTAU 

THIS  SUBROUTINE  COMPUTES  THE  ACCELERATION 
OF  THE  CELL  DUE  TO  THE  STRESSES 
IN  THE  RADIAL  DIRECTION. 

ACTING  ON  THIS  CELL  15  OF  THEM). 

STORE  THE  RADIAL  COMPONENT  OF  THE 
ACCELERATION  IN  XI. 

WS-TAUI I I ' / P I DY 
I  FI  GAM) 1*2 » I 

1  WSD-1. 

WSE= 1. 

GO  TO  3 

2  WSO  =  XUI) 

HSE=XI 1 1-1 ) 

3  WSA=DY(  JN)*(  SNR*WSD-SNL*'..3E) 

I  F ( GAM) 6.7.6 

6  WSF= I o 
GO  TO  4 

7  WSF=2® 

4  WSB=WS/WSF*ISTT~STB) 

WSC=DXI  1 1 ) Y { JN) *DDVK 
Xl=PIOY/AMXIK)*WSF*(  WSA-.  WSB-WSC) 

RETURN 
END 


IBFTC  DELTAV  LI  ST  .DECK , REF 
SUBROUTINE  DELTAV 
AXIAL  COMPONENT 

THIS  SUBROUTINE  COMPUTES  THE  ACCELERATION  OF 
THE  CELL  DUE  TO  THE  STRESSES  ACTING 
ON  THIS  CELL  (4  OF  THEM).  STORE  THE 
/  .IAL  COMPONENT  IN  X2. 

IF  (.GAM)  1.2.1 

1  WSB-1 . 

WSD=I • 

HSF=1. 

GO  TO  3 

2  WSB=2. 

WSD=XC  1 1 J 
WSF=X ( 1 1  — 1 ) 

3  WS=( SNT-SNB) /Wf  «TAUl II J/PIDY 
WSA=DY< JN «*(STR*WSD-STL*WSF) 
X2=PIDY/AMX(KJ*WSB*( WS+WSAJ 
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RETURN 

END 


$ IBFfC  ECALC  LI  ST, DECK. REF 

SUBROUTINE  ECALC 

THIS  ROUTINE  HILL  CALCULATE  THE  CHANGE 
IN  SPECIFIC  INTERNAL  ENERGY  DUE  TO  THE 
WORK  DONE  BY  THESE  STRESSES. 

STORE  IT  IN  (VI SCI. 

WSD=TAU( IiI*((U(K)+U(KA)J/2.*STT+ 
l(V{K)+V(KAi)/2.*SNT) 

WSE=TAU( I I )*l IU(KJ*U(KB1 1/2. *STB«- 
HV(K1*V(KB1  )/2.*SNBl 
IF ( GAR) 1,2,1 

1  WSF=l. 

MSG-1. 

WSH=1« 

GO  TO  3 

2  WSF=2. 

WSG=X (II 
MSH=X( 1 1— II 

3  HS=WSF*PIDY*DY( JNI 
WSA=MS*HSG*( (U(KR1*U(KJ 1/2.*SNR«- 

I(V(KR1+V(K1)/2#*STRI 
WSB=HS*WSH*( (U(K)+U(KLI I/2.^SNL+ 
i( VIK1+VIKLI 1/2.*STL1 
MSA=( WSA-WSB*MSD-WSE l/AHX(KJ*OT 
WSE=X1*QT 
MSD=X2^DT 

WSB=WSE*«  U(  K 1 +WSE/2. i+WSO* (  V I  K I  +MSD/2. 1 
VISC=WSA— MSB 
RETURN 
END 
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